home *** CD-ROM | disk | FTP | other *** search
/ Language/OS - Multiplatform Resource Library / LANGUAGE OS.iso / lisp / elk-2_0.lha / elk-2.0 / contrib / zelk / src-zelk / fvector.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-11-13  |  66.1 KB  |  2,880 lines
  1. /* fvector.c zilla 9sep91 - elk foreign vector(ized) operations
  2.  *
  3.  * NOTE this file must be compiled with ANSI C - the ARITHCASE macro
  4.  * needs ANSI macros.
  5.  *
  6.  * vector meaning vectorized, not scheme vector.
  7.  * c.f G.Blelloch, Vector Models for Data Parallel Computing, MIT Press
  8.  *
  9.  * todo:
  10.  *
  11.  *-(v-extract arr start len) extracts a subsequence.
  12.  * this can be accomplished using index,+,[], but add if for efficiency.
  13.  * This is done easily with gather, stride 1!
  14.  *-generalized matrix multiplier/stepper which can vectorize
  15.  * a convolution or matrix multiplication-like operation.
  16.  *
  17.  * issues:
  18.  * Should v-*,etc. also take scalars, or do we need a compiler?
  19.  * See parlet.e for a discussion. In brief,
  20.  * 1. making all functions accept any combination of scalar/vector
  21.  * requires ugly programming--see ARITHOP--currently arith ops allow this
  22.  * but comparisons do not.
  23.  * 2. This dynamic typing is easy for an interpreter but hard for a compiler--
  24.  * It can be difficult for a compiler to tell whether a variable currently 
  25.  * contains a vector.  In the future we might be compiling into C.
  26.  * 
  27.  * big problem is garbage collection - discarding farrays will
  28.  * quickly cause garbage collection, though little will need to be
  29.  * collected except the recently allocated farrays.
  30.  * How to restrict gc to the farrays??
  31.  * 1. A generational GC would work well!
  32.  * 2. implement a minicompiler for vector expressions,
  33.  *    have this compiler automatically manage the vector storage.
  34.  * 3. when make_vector is called, if it cannot Get_Bytes_NoGC(),
  35.  *    have it somehow gc recently allocated vectors.  vectors
  36.  *    could be on separate heap, but i dont see any way of
  37.  *    establishing gcness without a full gc, other than reference
  38.  *    counting.
  39.  * 4. (Viewpoint).  although gc may be triggered after a very few
  40.  *    vector ops, these ops are the equivalent of many scalar ops.
  41.  *    dont worry about it...
  42.  * 5. if farray data is allocated with malloc instead of Get_Bytes,
  43.  *    the problem changes: vector data will not trigger GC, but
  44.  *    also will not be freed until a GC is triggered by other means.
  45.  *
  46.  * Dataparallel primitives from chatterjee/belloch/fisher,
  47.  * "Size&access inference for data-parallel programs",
  48.  * Sigplan91 conf on prog lang design & implementation:
  49.  * plus-scan +\ out[0]=0; out[i]=out[i-1]+in[i-1] 1..size
  50.  * dist  out[i] = in1 0..in2
  51.  * distv out[i] = in1 0..in2_siz
  52.  * permute perm out[in2[i]] = in1[i] 0..in1_siz
  53.  * bpermute     out[i] = in1[in2[i]]
  54.  * dpermute     out[i]=in3[i]; out[in2[i]] = in1[i]
  55.  *
  56.     This file is Copyright (C) 1991 John Lewis
  57.  
  58.     This file is free software; you can redistribute it and/or modify
  59.     it under the terms of the GNU General Public License as published by
  60.     the Free Software Foundation.
  61.  
  62.     This program is distributed in the hope that it will be useful,
  63.     but WITHOUT ANY WARRANTY; without even the implied warranty of
  64.     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  65.     GNU General Public License for more details.
  66.  
  67.     You should have received a copy of the GNU General Public License
  68.     along with this program; if not, write to the Free Software
  69.     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  70.  
  71.  ****NOTE THE ELK COPYING GC: ALL Object REFERENCES MUST BE GC_LINKED
  72.  ****ACROSS CALLS WHICH MAY ALLOCATE STORAGE.  ALL C VARIABLES WHICH 
  73.  ****ARE ASSIGNED FROM THE ADDRESS OF AN OBJECT MUST BE REASSIGNED
  74.  ****AFTER A GC.
  75.  *
  76.  * modified
  77.  * 12oct        use random() not rnd() for release version
  78.  * 26sep        mapcount
  79.  * 12jul        bugfix in elevated_array: (v-* 3. (v-index 2)) was
  80.  *              incorrect- took its type from the array.
  81.  * 7jul        various operations changed to preserve shape
  82.  * 8jun         arith ops overloaded with scalars
  83.  * 4jun         v-truncate
  84.  * 24may        bugfix in v-shape, v-array,-ref,-set!
  85.  * 18may        shift,head,tail
  86.  * 17may        atan2, rotate
  87.  * 14may        v-append now varargs; v-reverse
  88.  * 11may        CORRECTED GC; arith-op loop unroll; shape[0] becomes
  89.  *              minor; v-[] handles matrix.
  90.  * 13apr        more good stuff. select
  91.  * 6apr         shape,reshape,more bool ops,compress,reference,
  92.  *              scans,reductions
  93.  * 5mar         distribute can take farray to specify size
  94.  * 31jan        added scatter
  95.  * 30jan        fixed gather, moved v-truncate -> farray-int,
  96.  *              remove obsolete v-aref
  97.  * 20jan        added comparison ops, not
  98.  * 9jan         added gather - semiduplicate of aref
  99.  * 12dec        moved primdef from here to zilla.c
  100.  * 2oct         added v-aref
  101.  * 11sep        gc-protect, add truncate,fmodulo,functions
  102.  * 10sep        bugfix in div,sub: not commutative
  103.  */
  104.  
  105. /*(DOCINIT
  106.   (MANPAGE
  107.     "/usr/local/pub/man/man3/elkvector.3"
  108.     "elkvector"
  109.     "3"
  110.     "scheme vector operations"))
  111. */
  112.  
  113. /*(SECTION
  114.   "INTRODUCTION"
  115.   "These functions perform a subset of APL-like vector operations"
  116.   "on Elk foreign arrays."
  117.   "Most functions operate on integer and float but not character farrays;"
  118.   "convert character farrays using e.g. {\kwd farray-int} before use."
  119.   "Scheme vector operations are currently restricted to rank <= 2,"
  120.   "and not all operations work for rank 2."
  121.   "Boolean vectors are represented as vectors of integers;"
  122.   "`word parallelism' is not used yet."
  123.   "The vector package extends foreign arrays to several dimensions."
  124.   "Arrays are stored in column-minor order."
  125.   ""
  126.   "References:"
  127.   "S.Kamin, Programming Languages, Addison-Wesley 1990, Ch.3;"
  128.   "G.Blelloch, Vector Models for Data-Parallel Computing, MIT, 1990.")
  129. */
  130.  
  131. /*(SECTION "FUNCTIONS")
  132. */
  133.  
  134. #include <theusual.h>
  135. #include <scheme.h>
  136. #include <zelk.h>
  137. #include <constants.h>
  138. #include <assert.h>
  139.  
  140. #if ELKVECTOR
  141.  
  142. /* type of boolean vector.
  143.  * chose int rather that char because ints may be faster to manipulate,
  144.  * and because more v-ops work on ints than on chars.
  145.  */
  146. #define Vbool int4
  147. #define T_Vbool T_Fixnum
  148.  
  149. #define Sym_integer     Intern("integer")
  150. #define Sym_real        Intern("real")
  151. #define Sym_string      Intern("string")
  152.  
  153. /*forward*/ static int type_elevate Zproto((Object,Object));
  154.  
  155. /* assert that A,B are both vectors of the same shape. return length */
  156. /* helper */
  157. int4 v_conform Zproto((Object,Object));
  158. int4 v_conform(A,B)
  159.   Object A,B;
  160. {
  161.   bool bit1,bit2;
  162.   int i;
  163.   Farray *a=FARRAY(A),*b=FARRAY(B);
  164.  
  165.   Check_Type(A,T_Farray);
  166.   Check_Type(B,T_Farray);
  167.  
  168.   bit1 = ((a->len != b->len) || (a->ndim != b->ndim));
  169.  
  170.   for( i=0, bit2=FALSE; i < FARRAY_MAXDIM; i++ ) {
  171.     if (a->shape[i] != b->shape[i]) {
  172.       bit2 = TRUE;
  173.       break;
  174.     }
  175.   }
  176.  
  177.   if (bit1 || bit2) 
  178.       Primitive_Error("vector shape mismatch");
  179.   return( a->len );
  180. } /*v_conform*/
  181.  
  182.  
  183. /* for mixed array/scalar, mixed int/float/byte operations--
  184.    assume A,B, or both are arrays.
  185.    Return new array of similar shape and elevated type.
  186.    Error if neither is an array, or if shapes are different. */
  187.  
  188. static Object v_elevated_array(Object A,Object B)
  189. {
  190.   bool aisarray = (TYPE(A)==T_Farray);
  191.   bool bisarray = (TYPE(B)==T_Farray);
  192.   Object C;
  193.   Ztrace(("fvector v_elevated_array %d %d\n",aisarray,bisarray));
  194.  
  195.   if (aisarray && bisarray) {
  196.     int len = v_conform(A,B);
  197.     C = farray_make(type_elevate(A,B),len);
  198.     farray_copyshape(A,C);
  199.   }
  200.   else if (aisarray) {
  201.     C = farray_make(type_elevate(A,B),FARRAY(A)->len);
  202.     farray_copyshape(A,C);
  203.   }
  204.   else if (bisarray) {
  205.     C = farray_make(type_elevate(A,B),FARRAY(B)->len);
  206.     farray_copyshape(B,C);
  207.   }
  208.   else Panic("v_elevated_array: neither!");
  209.  
  210.   return C;
  211. } /*v_elevated_array*/
  212.  
  213.  
  214.  
  215. /* return unraveled length specified by shape. for error checking */
  216. static int4 v_shapelen Zproto((int4 [],int4));
  217. static int4
  218. v_shapelen(shape,ndim)
  219.   int4 shape[],ndim;
  220. {
  221.   register int i;
  222.   register int4 len = 1;
  223.  
  224.   for( i=0; i < ndim; i++ )  len *= shape[i];
  225.  
  226.   return len;
  227. } /*shapelen*/
  228.  
  229.  
  230. /*%%%%%%%%%%%%%%%% multidimensional arrays %%%%%%%%%%%%%%%%*/
  231. /*(SECTION "Array Primitives"
  232. )*/
  233.  
  234. /*(DOCENTRY
  235.   (USAGE "(v-array type dim)"
  236.   "Dim can be an integer, resulting in a one-dimensional array,"
  237.   "or a foreign array of integers, resulting in a multidimensional array."
  238.   "(v-array 'real (% 2 3)) returns an array of 2 rows by 3 columns."
  239.   ))
  240. */
  241.  
  242.  
  243. #define VARRAY    Pvarray, "v-array", 2,2,EVAL,
  244.  
  245. Object Pvarray(Type,Len)
  246.   Object Type,Len;
  247. {
  248.   int type,ltype;
  249.   Object F;
  250.   Farray *f;
  251.   Error_Tag = "v-array";
  252.  
  253.   if (Type == Sym_real)                type = T_Flonum;
  254.   else if (Type == Sym_integer)        type = T_Fixnum;
  255.   else if (Type == Sym_string)         type = T_String;
  256.   else Primitive_Error("bad type");
  257.  
  258.   ltype = TYPE(Len);
  259.   if ((ltype == T_Fixnum) || (ltype == T_Bignum)) {
  260.     F = farray_make(type,Get_Integer(Len));
  261.     f = FARRAY(F);
  262.   }
  263.  
  264.   else if (ltype == T_Farray) {
  265.     Farray *l = FARRAY(Len);
  266.     register int i,len;
  267.     GC_Node;
  268.     if (l->type != T_Fixnum) Primitive_Error("bad array dimension");
  269.     len = 1;
  270.     for( i=0; i < l->len; i++ )  len *= ((int4 *)l->data)[i];
  271.  
  272.     GC_Link(Len);
  273.     F = farray_make(type,len);
  274.     GC_Unlink;
  275.  
  276.     f = FARRAY(F);
  277.     l = FARRAY(Len);
  278.     f->ndim = l->len;
  279.  
  280.     for( i=0; i < l->len; i++ ) {
  281.       /* note reversal of direction for storage */
  282.       /* in C, dimension[0] is column, but in scheme col is */
  283.       /* last=highest dimension */
  284.       f->shape[i] = ((int4 *)l->data)[(l->len-i)-1];
  285.     }
  286.   }
  287.  
  288.   else Primitive_Error("bad array dimension");
  289.  
  290.   /* Initialize arrays created with this (farray-make) primitive */
  291.   Zbzero((char *)f->data,((f->type)==T_String ? (f->len) : ((f->len)*4)));
  292.  
  293.   return F;
  294. } /*varray*/
  295.  
  296.  
  297.  
  298. /*(DOCENTRY
  299.   (USAGE "(v-array-set! v idx val)"
  300.   "Multidimensional accessing occurs when idx is a foreign array of integers."
  301.   "The indices are in the `normal' order: "
  302.   "the last element of a 2 (rows) by 3 array is (\% 1 2)."
  303.   "Since the index is itself an array, and due to the normal read/print order"
  304.   "for arrays, this means that the most major axis of the index is in [0]"
  305.   "and the minor axis is in the highest slot."
  306.   "Idx may also be an integer, in which case this is the same as farray-set!"
  307.   ))
  308. */
  309.  
  310.   /* helper- get a 1-D array offset from possibly multi-dimensional index */
  311.   static int vmdim_index Zproto((Farray *,Object));
  312.   static int vmdim_index(a,pidx)
  313.     Farray *a;
  314.     Object pidx;
  315.   {
  316.     int idx,ltype;
  317.  
  318.     ltype = (TYPE(pidx));
  319.     if ((ltype == T_Fixnum) || (ltype == T_Bignum)) {
  320.       idx = Get_Integer(pidx);
  321.     }
  322.  
  323.     else if (ltype == T_Farray) {
  324.       Farray *aidx = FARRAY(pidx);
  325.       int i,ii,axisstep;
  326.       if ((aidx->len != a->ndim) ||
  327.           (aidx->type != T_Fixnum)) Primitive_Error("bad index");
  328.       idx = 0;
  329.       axisstep = 1;  /* e.g. z*(r*c) + y*c + x */
  330.       
  331.       for( i=0; i < a->ndim; i++ ) {
  332.         ii = ((int *)aidx->data)[(a->ndim - i)-1];  /* reverse order! */
  333.         if ((ii < 0) || (ii >= a->shape[i]))
  334.           Primitive_Error("index out of range");
  335.         if (i > 0) axisstep *= a->shape[i-1];
  336.         idx += ii * axisstep;
  337.       }
  338.     }
  339.     else Primitive_Error("bad index type"); 
  340.  
  341.     if ((idx < 0) || (idx >= a->len)) Primitive_Error("index out of array");
  342.  
  343.     return idx;
  344.   } /*vmdim_index*/
  345.  
  346. #define VARRAYSET    Pvarray_set, "v-array-set!", 3,3,EVAL,
  347.  
  348. Object Pvarray_set(f,pidx,pobj)
  349.   Object f,pidx,pobj;
  350. {
  351.   Farray *a;
  352.   int4 idx;
  353.   long *L; float *F; unsigned char *C;
  354.  
  355.   Error_Tag = "v-array-set!";
  356.   Check_Type(f,T_Farray);
  357.  
  358.   a = FARRAY(f);
  359.   C = (unsigned char *)a->data;
  360.   F = (float *)a->data;
  361.   L = (long *)a->data;
  362.   
  363.   idx = vmdim_index(a,pidx);
  364.  
  365.   switch(a->type) {
  366.   case T_Fixnum:
  367.     L[idx] = Get_Integer(pobj); 
  368.     break;
  369.   case T_Flonum: 
  370.     if (TYPE(pobj) != T_Flonum)  Primitive_Error("bad type");
  371.     F[idx] = (double)FLONUM(pobj)->val;
  372.     break;
  373.   case T_String:
  374.     /*    if (TYPE(pobj) != T_Character)  Primitive_Error("bad type"); 
  375.           C[idx] = (char)CHAR(pobj); */
  376.     C[idx] = (unsigned char)Get_Integer(pobj);
  377.     break;
  378.   default: Panic("varray_set");
  379.   } /*switch a->type*/
  380.  
  381.   return pobj;
  382. } /*P_set*/
  383.  
  384.  
  385. /*(DOCENTRY
  386.   (USAGE "(v-array-ref v idx)"
  387.   "Multidimensional accessing occurs when idx is a foreign array of integers."
  388.   "The indices are in the `normal' order: "
  389.   "the last element of a 2 (rows) by 3 array is (% 1 2)."
  390.   "Since the index is itself an array, and due to the normal read/print order"
  391.   "for arrays, this means that the most major axis of the index is in [0]"
  392.   "and the minor axis is in the highest slot."
  393.   "Idx may also be an integer, in which case this is the same as farray-ref"
  394.   ))
  395. */
  396.  
  397. #define VARRAYREF    Pvarray_ref, "v-array-ref", 2,2,EVAL,
  398.  
  399. Object Pvarray_ref(f,pidx)
  400.   Object f,pidx;
  401. {
  402.   int4 idx;
  403.   Farray *a;
  404.   long *L; float *F; unsigned char *C;
  405.   Object val;
  406.   Error_Tag = "v-array-ref";
  407.  
  408.   Check_Type(f,T_Farray);
  409.  
  410.   a = FARRAY(f);
  411.   C = (unsigned char *)a->data;
  412.   F = (float *)a->data;
  413.   L = (long *)a->data;
  414.  
  415.   idx = vmdim_index(a,pidx);
  416.  
  417.   switch(a->type) {
  418.   case T_Fixnum:
  419.     val = Make_Integer((int4)L[idx]);
  420.     break;
  421.   case T_Flonum:
  422.     val = Make_Reduced_Flonum(F[idx]);
  423.     break;
  424.   case T_String:
  425. /*  val = Make_Char(C[idx]); */
  426.     val = Make_Integer((int4)C[idx]);
  427.     break;
  428.   default: Panic("farray_ref");
  429.   }
  430.  
  431.   return val;
  432. } /*varray_ref*/
  433.  
  434.  
  435. /*%%%%%%%%%%%%%%%% shaping %%%%%%%%%%%%%%%%*/
  436. /*(SECTION "Shape Functions")
  437. */
  438.  
  439. /*(DOCENTRY
  440.   (USAGE "(v-shape v) => returns the length or shape of v in an int farray."
  441.   "For example, the shape of a 2 (rows) by 3 matrix is (% 2 3)."
  442.   "BEWARE: "
  443.   "Since arrays are printed from slot 0 up, the major axis is actually"
  444.   "stored in slot 0 of the result, and the minor axis is in the highest slot."
  445.   "V-shape can thus be passed to v-array directly."
  446.   ))
  447. */
  448.  
  449. #define VSHAPE    Pvshape, "v-shape", 1,1,EVAL,
  450. Object 
  451. Pvshape(A)
  452.   Object A;
  453. {
  454.   register int i;
  455.   Object B;
  456.   Farray *a;
  457.   register int4 *ib;
  458.   GC_Node;
  459.   Error_Tag = "v-shape";
  460.  
  461.   Check_Type(A,T_Farray);
  462.  
  463.   GC_Link(A);
  464.   B = farray_make(T_Fixnum,FARRAY(A)->ndim);
  465.   GC_Unlink;
  466.  
  467.   a = FARRAY(A);
  468.  
  469.   ib = (int4 *)FARRAY(B)->data;
  470.   for( i=0; i < a->ndim; i++ )
  471.     *ib++ = a->shape[(a->ndim-i)-1];
  472.  
  473.   return B;
  474. } /*Pvshape*/
  475.  
  476.  
  477. /*(DOCENTRY
  478.   (USAGE "(v-ravel v) => returns an unraveled vector,"
  479.          "ie, vector containing concatenation of matrix rows or planes"))
  480. */
  481. #define VRAVEL    Pvravel, "v-ravel", 1,1,EVAL,
  482.  
  483. static Object
  484. Pvravel(A)
  485.   Object A;   
  486. {
  487.   Object B;
  488.   Farray *b;
  489.  
  490.   B = P_farray_copy(A);
  491.   b = FARRAY(B);
  492.   b->ndim = 1;
  493.   b->shape[0] = b->len;
  494.  
  495.   return B;
  496. } /*ravel*/
  497.  
  498.  
  499. /*(DOCENTRY
  500.   (USAGE "(reshape array shape[ndim]) --  Give a new shape to an array."
  501.          "Sometimes called restruct."
  502.          "In Kamin this is defined to fill the specified shape, "
  503.          "wrapping if needed."
  504.          "The specified shape and source length conform,"
  505.          "whereas in Kamin reshape wraps if needed."))
  506. */
  507. #define VRESHAPE    Pvreshape, "v-reshape", 2,2,EVAL,
  508. static Object Pvreshape(A,Shape)
  509.   Object A,Shape;
  510. {
  511.   int i;
  512.   Object B;
  513.   Farray *s,*b;
  514.   register int4 *is;
  515.   GC_Node2;
  516.   Error_Tag = "v-reshape";
  517.  
  518.   Check_Type(A,T_Farray);
  519.   Check_Type(Shape,T_Farray);
  520.  
  521.   s = FARRAY(Shape);
  522.   if (s->type != T_Fixnum)
  523.     Primitive_Error("shape must be array of fixnums");
  524.  
  525.   if ((s->len < 0) || (s->len > FARRAY_MAXDIM))
  526.     Primitive_Error("bad # of dimensions");
  527.  
  528.   if (v_shapelen((int4 *)s->data,s->len) != FARRAY(A)->len)
  529.     Primitive_Error("shape does not fit array");
  530.  
  531.   GC_Link2(A,Shape);  /* protecting A not necessary here */
  532.   B = P_farray_copy(A);
  533.   GC_Unlink;
  534.  
  535.   b = FARRAY(B);
  536.   s = FARRAY(Shape);    /* reassign - Shape may have moved during gc */
  537.   b->ndim = s->len;
  538.  
  539.   is = (int4 *)s->data;
  540.   for( i=b->ndim-1; i >= 0; i-- )
  541.     b->shape[i] = *is++;
  542.  
  543.   return B;
  544. } /*Pvreshape*/
  545.  
  546.  
  547.  
  548. /*(DOCENTRY
  549.   (USAGE "(v-transpose v) => Return transpose of matrix."))
  550. */
  551.  
  552. #define VTRANSPOSE    Pvtranspose, "v-transpose", 1,1,EVAL,
  553.  
  554. /* nth column becomes nth row */
  555. static Object Pvtranspose(A)
  556.   Object A;
  557. {
  558.   Object B;
  559.   Farray *a,*b;
  560.   register int nrows,ncols;
  561.   GC_Node;
  562.   Error_Tag = "v-transpose";
  563.  
  564.   Check_Type(A,T_Farray);
  565.   if (FARRAY(A)->ndim != 2) Primitive_Error("arg is not a matrix");
  566.  
  567.   a = FARRAY(A);
  568.   GC_Link(A);
  569.   B = farray_make(a->type,a->len);
  570.   GC_Unlink;
  571.  
  572.   a = FARRAY(A);        /* reassign - A may have moved during GC */
  573.   b = FARRAY(B);
  574.   b->ndim = 2;
  575.   b->shape[0] = nrows = a->shape[1];
  576.   b->shape[1] = ncols = a->shape[0];
  577.   Ztrace(("v-transpose [%d,%d]\n",nrows,ncols));
  578.  
  579. #define XPOSE(type) \
  580.   { \
  581.     register type *ap = (type *)a->data;\
  582.     register type *bp = (type *)b->data;\
  583.     register int r,c;\
  584.     for( r=0; r < nrows; r++ ) {\
  585.       for( c=0; c < ncols; c++ ) {\
  586.         bp[c*nrows+r] = ap[r*ncols+c];\
  587.       }\
  588.     }\
  589.   }
  590.  
  591.   if (a->type == T_Flonum) 
  592.     XPOSE(float)
  593.  
  594.   else if (a->type == T_Fixnum)
  595.     XPOSE(int4)
  596.  
  597.   else if (a->type == T_String)
  598.     XPOSE(unsigned char)
  599.  
  600.   else Panic("transpose-datatype?");
  601.   
  602.   return B;
  603. } /*transpose*/
  604.  
  605.  
  606.  
  607. /*(DOCENTRY
  608.   (USAGE "(v-compress v boolvector) => Return vector of elements of v"
  609.          "corresponding to ones in boolvector."
  610.          "Called compress in APL, pack in Blelloch."))
  611. */
  612.  
  613. /* helper to compress */
  614. static Object vcompress_matrix Zproto((Object,Object));
  615. static Object vcompress_matrix(A,Bitvec)
  616.   Object A,Bitvec;
  617. {
  618.   Primitive_Error("not implemented for non-vector");
  619. }
  620.  
  621.  
  622. /* helper to compress */
  623. static Object vcompress_vector Zproto((Object,Object));
  624. static Object vcompress_vector(A,Bitvec)
  625.   Object A,Bitvec;
  626. {
  627.   Farray *a,*bv;
  628.   Object B;
  629.   Farray *b;
  630.   register int4 *ib;
  631.   register int i,bvlen,sum;
  632.   GC_Node2;
  633.  
  634. /* count up the number of set bits = length of compressed vector*/
  635.   bv = FARRAY(Bitvec);
  636.   a = FARRAY(A);
  637.   ib = (Vbool *)bv->data;
  638.   bvlen = bv->len;
  639.   sum = 0;
  640.   for( i=0; i < bvlen; i++ )  if (*ib++) sum++;
  641.  
  642. /* copy to new array */
  643.   GC_Link2(A,Bitvec);
  644.   B = farray_make(a->type,sum);
  645.   GC_Unlink;
  646.  
  647.   b = FARRAY(B);        
  648.   a = FARRAY(A);        /* reassign - A,Bitvec may have moved */
  649.   bv = FARRAY(Bitvec);
  650.  
  651. # define COMPRESSLOOP(type) {\
  652.     type *bp,*ap;\
  653.     bp = (type *)b->data;\
  654.     ap = (type *)a->data;\
  655.     ib = (Vbool *)bv->data;\
  656.     for( i=0; i < bvlen; i++ ) {\
  657.       if (*ib++) *bp++ = *ap;\
  658.       ap++;\
  659.     }\
  660.   }
  661.  
  662.   switch(a->type) {
  663.   case T_Flonum:
  664.     COMPRESSLOOP(float);
  665.     break;
  666.  
  667.   case T_Fixnum:
  668.     COMPRESSLOOP(int4);
  669.     break;
  670.  
  671.   case T_String:
  672.     COMPRESSLOOP(Zbyte);
  673.     break;
  674.  
  675.   default: Panic("v-compress");
  676.     break;
  677.   }
  678. # undef COMPRESSLOOP 
  679.   return B;
  680. } /*vcompress_vector*/
  681.  
  682.  
  683.  
  684. #define VCOMPRESS    Pvcompress, "v-compress", 2,2,EVAL,
  685. Object 
  686. Pvcompress(A,Bitvec)
  687.   Object A,Bitvec;
  688. {
  689.   Object B;
  690.   Farray *a,*bv;
  691.   Error_Tag = "v-compress";
  692.  
  693.   Check_Type(A,T_Farray);
  694.   Check_Type(Bitvec,T_Farray);
  695.  
  696.   a = FARRAY(A);
  697.   bv = FARRAY(Bitvec);
  698.   if (bv->type != T_Vbool)
  699.     Primitive_Error("arg2 must be boolean vector");
  700.  
  701.   if (a->ndim != 1) {
  702.     if (a->shape[1] != bv->len)
  703.       Primitive_Error("length of arg2 must match major axis of arg1");
  704.     B = vcompress_matrix(A,Bitvec);
  705.   }
  706.   else {
  707.     if (a->len != bv->len)
  708.       Primitive_Error("length of arg2 must match length of arg1");
  709.     B = vcompress_vector(A,Bitvec);
  710.   }
  711.  
  712.   return B;
  713. } /*compress*/
  714.  
  715.  
  716.  
  717. /*(MANENTRY
  718.   "(v-append v1 v2) => Concatenate v1,v2")
  719. */
  720.  
  721.  
  722. #define VAPPEND    Pvappend, "v-append", 0,MANY,VARARGS,
  723. Object 
  724. Pvappend(ac,av)
  725.   int ac;     
  726.   Object *av;
  727. {
  728.   Object C;
  729.   int ic,clen;
  730.   Error_Tag = "v-append";
  731.  
  732.   clen = 0;
  733.   for( ic=0; ic < ac; ic++) {
  734.     Check_Type(av[ic],T_Farray);
  735.     if (FARRAY(av[ic])->type != FARRAY(av[0])->type)
  736.       Primitive_Error("types must match");
  737.     clen += FARRAY(av[ic])->len;
  738.   }
  739.  
  740.   C = farray_make(FARRAY(av[0])->type,clen);
  741.  
  742. # define COPYTYPE(type) {\
  743.     register type *x = (type *)FARRAY(C)->data;\
  744.         \
  745.     for( ic=0; ic < ac; ic++ ) {\
  746.        register type *y;\
  747.        register int4 i,len;\
  748.         \
  749.        y = (type *)FARRAY(av[ic])->data;\
  750.        len = FARRAY(av[ic])->len;\
  751.        for( i=0; i < len; i++ ) *x++ = *y++;\
  752.      }\
  753.   }\
  754.     
  755.   switch(FARRAY(av[0])->type) {
  756.   case T_String: COPYTYPE(unsigned char) break; 
  757.   case T_Fixnum: COPYTYPE(int) break; 
  758.   case T_Flonum: COPYTYPE(float) break;
  759.   default: Panic("v-append");
  760.   }
  761. # undef COPYTYPE
  762.     
  763.   return C;
  764. } /*append*/
  765.  
  766.  
  767.  
  768. /*(MANENTRY
  769.   "(v-select b v1 v2) => select between v1,v2."
  770.   "Copy v1[i] if b[i]=1, else v2[i].")
  771. */
  772.  
  773. #define VSELECT    Pvselect, "v-select", 3,3,EVAL,
  774. Object 
  775. Pvselect(B,V1,V2)
  776.   Object B,V1,V2;
  777. {
  778.   Object V3;
  779.   Farray *b;
  780.   GC_Node3;
  781.   register int i,len;
  782.   Error_Tag = "v-select";
  783.  
  784.   len = v_conform(V1,V2);
  785. /*  Check_Type(V1,T_Farray); not necessary - v_conform checks
  786.    Check_Type(V2,T_Farray); */
  787.  
  788.   b = FARRAY(B);
  789.   if (b->type != T_Vbool)
  790.     Primitive_Error("bool vector required");
  791.   if (b->len != len)
  792.     Primitive_Error("length mismatch");
  793.  
  794.   GC_Link3(B,V1,V2);
  795.   V3 = farray_make(FARRAY(V1)->type,FARRAY(V1)->len);
  796.   GC_Unlink;
  797.  
  798.   b = FARRAY(B);
  799.  
  800. # define SELECT(vtype) \
  801.   {\
  802.     vtype *vp1 = (vtype *)FARRAY(V1)->data;\
  803.     vtype *vp2 = (vtype *)FARRAY(V2)->data;\
  804.     vtype *vp3 = (vtype *)FARRAY(V3)->data;\
  805.     Vbool *vpb = (Vbool *)b->data;\
  806.     for( i=0; i < len; i++ ) {\
  807.       *vp3 = *vpb ? *vp1 : *vp2;\
  808.       vp3++; vp1++; vp2++; vpb++;\
  809.     }\
  810.   }
  811.  
  812.   switch(FARRAY(V1)->type) {
  813.   case T_String: SELECT(unsigned char) break;
  814.   case T_Fixnum: SELECT(int4) break;
  815.   case T_Flonum: SELECT(float) break;
  816.   default: Panic("v-select");
  817.   }
  818. # undef SELECT
  819.     
  820.   return V3;
  821. } /*select*/
  822.  
  823.  
  824. /*%%%%%%%%%%%%%%%% gather/scatter/subscripting %%%%%%%%%%%%%%%%*/
  825. /*(SECTION "Gather/Scatter/Subscripting")
  826. */
  827.  
  828.  
  829. /*(DOCENTRY
  830.   (USAGE "(v-mapcount v min max) -- vector map"
  831.          "Returns an array min..max which counts the number of times"
  832.          "a particular value appears in v. Values outside min..max"
  833.          "are discarded.  INTEGER ONLY!"
  834.          "Useful for histogram generation for example."
  835.          "Similar to map primitive in Paragon")
  836.   )
  837. */
  838.  
  839. #define VMAPCOUNT    Pvmapcount, "v-mapcount", 3,3,EVAL,
  840.  
  841. static Object Pvmapcount Zproto((Object,Object,Object));
  842. static Object Pvmapcount(V,Min,Max)
  843.   Object V,Min,Max;
  844. {
  845.   Object Rval;
  846.   register int4 *r,*v;
  847.   register int i,vlen,min,max;
  848.   int rlen;
  849.   GC_Node;
  850.   Error_Tag = "v-mapcount";
  851.  
  852.   min = Get_Integer(Min);
  853.   max = Get_Integer(Max);
  854.   rlen = 1 + (max - min);
  855.   if (rlen < 1) Primitive_Error("bad min,max");
  856.  
  857.   Check_Type(V,T_Farray);
  858.   if (FARRAY(V)->type != T_Fixnum) Primitive_Error("integer only");
  859.  
  860.   GC_Link(V);
  861.   Rval = farray_make(T_Fixnum,rlen);
  862.   GC_Unlink;
  863.  
  864.   r = (int4 *)FARRAY(Rval)->data;
  865.   Zbzero((char *)FARRAY(Rval)->data,sizeof(int4) * rlen);
  866.  
  867.   v = (int4 *)FARRAY(V)->data;
  868.   vlen = FARRAY(V)->len;
  869.  
  870.   for( i=0; i < vlen; i++ ) {
  871.     register int iv = v[i];
  872.     if ((iv < min) || (iv > max)) continue;
  873.     r[iv-min] ++;
  874.   }
  875.  
  876.   return Rval;
  877. } /*mapcount*/
  878.  
  879.  
  880. /*(DOCENTRY
  881.   (USAGE "(v-[] v p) -- subscripting"
  882.          "If p is scalar, this is equivalent to farray-ref."
  883.          "If p is a vector, the corresponding elements are"
  884.          "returned from v (the position elements need not be unique),"
  885.          "and the result has the length of p."
  886.          "If v is a matrix, the rows corresponding to p are returned. "
  887.          "Source--Kamin; c.f. permute in Blelloch p.62, which is different."))
  888. */
  889.  
  890. /* helper to reference */
  891. static Object vreference_matrix Zproto((Object,Object));
  892. static Object vreference_matrix(A,Ref)
  893.   Object A,Ref;
  894. {
  895.   Object B;
  896.   int blen;
  897.   Farray *a,*r,*b;
  898.   int ir;
  899.   int ncols,nrows;
  900.   GC_Node2;
  901.  
  902.   if (FARRAY(A)->ndim != 2) Panic("vref_mat");
  903.  
  904.   /* # of rows * # of columns */
  905.   blen = FARRAY(Ref)->len * FARRAY(A)->shape[0];
  906.  
  907.   GC_Link2(A,Ref);
  908.   B = farray_make(FARRAY(A)->type,blen);
  909.   GC_Unlink;
  910.  
  911.   r = FARRAY(Ref);
  912.   a = FARRAY(A);
  913.   b = FARRAY(B);
  914.  
  915.   ncols = b->shape[0] = a->shape[0];
  916.   nrows = b->shape[1] = r->len;
  917.   b->ndim = 2;
  918.  
  919. # define REFCOPY(typ_) {\
  920.     int rowsize = ncols * sizeof(typ_); \
  921.     typ_ *bp = (typ_ *)b->data; \
  922.     for( ir=0; ir < nrows; ir++ ) { \
  923.       typ_ *ap = ((typ_ *)a->data) + ((int4 *)r->data)[ir]*ncols; \
  924.       Zbcopy((Zunspec)ap,(Zunspec)bp,rowsize); \
  925.       bp += ncols; \
  926.     } \
  927.   }
  928.  
  929.   switch(a->type) {
  930.   case T_Fixnum: REFCOPY(int4); break;
  931.  
  932.   case T_Flonum: REFCOPY(float); break;
  933.  
  934.   case T_String: REFCOPY(unsigned char *); break;
  935.  
  936.   default: Panic("vref_mat(2)");
  937.   }
  938.     
  939.   return B;
  940. # undef REFCOPY
  941. } /*vreference_matrix*/
  942.  
  943.  
  944. /* helper to reference */
  945. static Object vreference_vector Zproto((Object,Object));
  946. static Object vreference_vector(A,Ref)
  947.   Object A,Ref;
  948. {
  949.   Farray *a,*r;
  950.   Object B;
  951.   Farray *b;
  952.   GC_Node2;
  953.  
  954.   GC_Link2(A,Ref);
  955.   /* NB type of a, length of r */
  956.   B = farray_make(FARRAY(A)->type,FARRAY(Ref)->len);
  957.   GC_Unlink;
  958.  
  959.   a = FARRAY(A);
  960.   r = FARRAY(Ref);
  961.   b = FARRAY(B);
  962.  
  963. # define REFERENCELOOP(type) {\
  964.     register int i;\
  965.     register int4 *rp;\
  966.     register int rlen = r->len;\
  967.     register int alen = a->len;\
  968.     type *bp,*ap;\
  969.     bp = (type *)b->data;\
  970.     ap = (type *)a->data;\
  971.     rp = (int4 *)r->data;\
  972.     for( i=0; i < rlen; i++ ) {\
  973.       if ((*rp < 0) || (*rp >= alen))\
  974.         Primitive_Error("reference out of vector");\
  975.       *bp++ = ap[*rp++];\
  976.     }\
  977.   }
  978.  
  979.   switch(a->type) {
  980.   case T_Flonum:
  981.     REFERENCELOOP(float);
  982.     break;
  983.  
  984.   case T_Fixnum:
  985.     REFERENCELOOP(int4);
  986.     break;
  987.  
  988.   case T_String:
  989.     REFERENCELOOP(Zbyte);
  990.     break;
  991.  
  992.   default: Panic("v-reference");
  993.     break;
  994.   }
  995. # undef REFERENCELOOP 
  996.   return B;
  997. } /*vreference_vector*/
  998.  
  999.  
  1000.  
  1001. #define VREFERENCE    Pvreference, "v-reference", 2,2,EVAL,
  1002. #define VREFERENCEb    Pvreference, "v-[]", 2,2,EVAL,
  1003. Object 
  1004. Pvreference(A,Ref)
  1005.   Object A,Ref;
  1006. {
  1007.   Object B;
  1008.   Farray *a,*r;
  1009.   Error_Tag = "v-reference";
  1010.  
  1011.   Check_Type(A,T_Farray);
  1012.  
  1013.   if (TYPE(Ref) != T_Farray)  /* simple scalar reference */
  1014.     return P_farray_ref(A,Ref);
  1015.  
  1016.   a = FARRAY(A);
  1017.   r = FARRAY(Ref);
  1018.   if (r->type != T_Fixnum)
  1019.     Primitive_Error("subscripts must be integer");
  1020.  
  1021.   if (a->ndim != 1) {
  1022.     B = vreference_matrix(A,Ref);
  1023.   }
  1024.   else {
  1025. /*  this restriction is not necessary:
  1026.     if (a->len < r->len)
  1027.       Primitive_Error("length of arg2 must be LE length of arg1");
  1028.  */
  1029.     B = vreference_vector(A,Ref);
  1030.   }
  1031.  
  1032.   return B;
  1033. } /*reference*/
  1034.  
  1035.  
  1036.  
  1037.  
  1038. /*(DOCENTRY
  1039.   (USAGE "(v-gather v offset stride len)"
  1040.          "Return a row or column of a multidimensional array,"
  1041.          "or similar operations."
  1042.          "   gather[j] = source[offset + j*stride];  j=0..len-1"))
  1043. */
  1044.  
  1045. /* "With a row-major 2D array, consider:"
  1046.    " Returning column c requires offset=c, stride=ncols, len=nrows."
  1047.    "  len could be determined automatically."
  1048.    " Returning row r requires offset=r*ncols, stride=1, len=ncols"
  1049.    "  len cannot be determined from offset and stride alone:"
  1050.    "  consider 3rd row (counting from 1) of a 4x4 array:"
  1051.    "  offset = 2*4=8, stride=1.  We need to know len(=ncols),"
  1052.    "  and we know only that the ncols divides offset evenly - (2? 4? 8?)."
  1053.    "  > Thus, len must be provided."
  1054. */
  1055.  
  1056. #define VGATHER    Pvgather, "v-gather", 4,4,EVAL,
  1057. Object Pvgather(F,Offset,Stride,Len)
  1058.   Object F,Offset,Stride,Len;
  1059. {
  1060.   Farray *fa,*ga;
  1061.   Object G;
  1062.   int stride,offset,len;
  1063.   int gi;
  1064.   GC_Node;
  1065.   Error_Tag = "v-gather";
  1066.  
  1067.   Check_Type(F,T_Farray);
  1068.  
  1069.   offset = Get_Integer(Offset);
  1070.   stride = Get_Integer(Stride);
  1071.   len = Get_Integer(Len);
  1072.  
  1073.   if (stride <= 0)
  1074.     Primitive_Error("bad stride");
  1075.   if ((offset < 0) || (offset >= FARRAY(F)->len))
  1076.     Primitive_Error("bad offset");
  1077.   if (len <= 0)
  1078.     Primitive_Error("bad len");
  1079.  
  1080.   if ((offset + (len-1)*stride) >= FARRAY(F)->len)
  1081.     Primitive_Error("offset+len*stride is beyond array");
  1082.  
  1083.   GC_Link(F);
  1084.   G = farray_make(FARRAY(F)->type,len);
  1085.   GC_Unlink;
  1086.  
  1087.   fa = FARRAY(F);
  1088.   ga = FARRAY(G);
  1089.  
  1090.   switch(fa->type) {
  1091.  
  1092.   case T_Flonum: {
  1093.     float *fd,*gd;
  1094.     fd = ((float *)fa->data) + offset;
  1095.     gd = (float *)ga->data;
  1096.     for( gi=0; gi < len; gi++ ) {
  1097.       *gd++ = *fd;
  1098.       Ztrace(("farray-gather [%d]->[%d]\n",
  1099.               fd-(float *)fa->data,(gd-1)-(float *)ga->data));
  1100.       fd += stride;
  1101.     }
  1102.     /* assert that last iteration did not walk off end of array */
  1103.     assert( (fd - stride) < (((float *)fa->data) + fa->len) );
  1104.     break;
  1105.   }
  1106.  
  1107.   case T_Fixnum: {
  1108.     int *fd,*gd;
  1109.     fd = ((int *)fa->data) + offset;
  1110.     gd = (int *)ga->data;
  1111.     for( gi=0; gi < len; gi++ ) {
  1112.       *gd++ = *fd;
  1113.       fd += stride;
  1114.     }
  1115.     /* assert that last iteration did not walk off end of array */
  1116.     assert( (fd - stride) < (((int *)fa->data) + fa->len) );
  1117.     break;
  1118.   }
  1119.  
  1120.   case T_String: {
  1121.     unsigned char *fd,*gd;
  1122.     fd = ((unsigned char *)fa->data) + offset;
  1123.     gd = (unsigned char *)ga->data;
  1124.     for( gi=0; gi < len; gi++ ) {
  1125.       *gd++ = *fd;
  1126.       fd += stride;
  1127.     }
  1128.     /* assert that last iteration did not walk off end of array */
  1129.     assert( (fd - stride) < (((unsigned char *)fa->data) + fa->len) );
  1130.     break;
  1131.   }
  1132.  
  1133.   } /*switch*/
  1134.  
  1135.   return G;
  1136. } /*gather*/
  1137.  
  1138.  
  1139.  
  1140. /*(DOCENTRY
  1141.   (USAGE "(v-scatter destarray offset stride len srcarray)"
  1142.          "Inverse of gather."
  1143.          "Insert a row or column into a multidimensional array,"
  1144.          "or similar operations."
  1145.          "   src[j] => dest[offset + j*stride];  j=0..len-1"
  1146.          "NOTE THIS IS A SIDE EFFECT--DATA IS INSERTED INTO EXISTING DEST"))
  1147. */
  1148.  
  1149. #define VSCATTER    Pvscatter, "v-scatter", 5,5,EVAL,
  1150. Object Pvscatter(Dest,Offset,Stride,Len, Src)
  1151.   Object Dest,Offset,Stride,Len,Src;
  1152. {
  1153.   Farray *da,*sa;
  1154.   int stride,offset,len;
  1155.   int i;
  1156.   Error_Tag = "v-scatter";
  1157.  
  1158.   Check_Type(Dest,T_Farray);
  1159.   Check_Type(Src,T_Farray);
  1160.  
  1161.   da = FARRAY(Dest);
  1162.   sa = FARRAY(Src);
  1163.   offset = Get_Integer(Offset);
  1164.   stride = Get_Integer(Stride);
  1165.   len = Get_Integer(Len);
  1166.  
  1167.   if (stride <= 0)
  1168.     Primitive_Error("bad stride");
  1169.   if ((offset < 0) || (offset >= da->len))
  1170.     Primitive_Error("bad offset");
  1171.   if (len <= 0)
  1172.     Primitive_Error("bad len");
  1173.  
  1174.   if ((offset + (len-1)*stride) >= da->len)
  1175.     Primitive_Error("offset+len*stride is beyond array");
  1176.  
  1177.   if (da->type != sa->type)
  1178.     Primitive_Error("array type mismatch");
  1179.  
  1180.   switch(da->type) {
  1181.  
  1182.   case T_Flonum: {
  1183.     float *dd,*sd;
  1184.     dd = ((float *)da->data) + offset;
  1185.     sd = (float *)sa->data;
  1186.     for( i=0; i < len; i++ ) {
  1187.       *dd = *sd++;
  1188.       Ztrace(("darray-scatter [%d]<-[%d]\n",
  1189.               dd-(float *)da->data,(sd-1)-(float *)sa->data));
  1190.       dd += stride;
  1191.     }
  1192.     /* assert that last iteration did not walk off end of array */
  1193.     assert( (dd - stride) < (((float *)da->data) + da->len) );
  1194.     break;
  1195.   }
  1196.  
  1197.   case T_Fixnum: {
  1198.     int *dd,*sd;
  1199.     dd = ((int *)da->data) + offset;
  1200.     sd = (int *)sa->data;
  1201.     for( i=0; i < len; i++ ) {
  1202.       *dd = *sd++;
  1203.       dd += stride;
  1204.     }
  1205.     /* assert that last iteration did not walk off end of array */
  1206.     assert( (dd - stride) < (((int *)da->data) + da->len) );
  1207.     break;
  1208.   }
  1209.  
  1210.   case T_String: {
  1211.     unsigned char *dd,*sd;
  1212.     dd = ((unsigned char *)da->data) + offset;
  1213.     sd = (unsigned char *)sa->data;
  1214.     for( i=0; i < len; i++ ) {
  1215.       *dd = *sd++;
  1216.       dd += stride;
  1217.     }
  1218.     /* assert that last iteration did not walk off end of array */
  1219.     assert( (dd - stride) < (((unsigned char *)da->data) + da->len) );
  1220.     break;
  1221.   }
  1222.  
  1223.   } /*switch*/
  1224.  
  1225.   return Dest;
  1226. } /*scatter*/
  1227.  
  1228.  
  1229.  
  1230. /*%%%%%%%%%%%%%%%% unary/monadic functions %%%%%%%%%%%%%%%%*/
  1231. /*(SECTION
  1232.   "Monadic Functions"
  1233.   "With the exceptions of v-abs and v-not, the result vectors are float.")
  1234. */
  1235.  
  1236.  
  1237. /* helper: apply function f to array A */
  1238. static Object vfuncall Zproto(( Object, double (*)(double) ));
  1239. static Object vfuncall(A,f)
  1240.   Object A;
  1241.   double (*f) Zproto((double));
  1242. {
  1243.   register int i,len;
  1244.   register float *ib;
  1245.   Object B;
  1246.   Farray *a;
  1247.   GC_Node;
  1248.  
  1249.   Check_Type(A,T_Farray);
  1250.  
  1251.   GC_Link(A);
  1252.   B = farray_make(T_Flonum,FARRAY(A)->len);
  1253.   GC_Unlink;
  1254.  
  1255.   a = FARRAY(A);
  1256.   len = a->len;
  1257.   ib = (float *)FARRAY(B)->data;
  1258.  
  1259.   switch(a->type) {
  1260.   case T_Fixnum: {
  1261.     register int4 *ia;
  1262.     ia = (int4 *)a->data;
  1263.     for( i=0; i < len; i++ ) *ib++ = (float)(*f)( (double)*ia++ );
  1264.     break;
  1265.   }
  1266.  
  1267.   case T_Flonum: {
  1268.     register float *ia;
  1269.     ia = (float *)a->data;
  1270.     for( i=0; i < len; i++ ) *ib++ = (*f)( *ia++ );
  1271.     break;
  1272.   }
  1273.  
  1274.   default:
  1275.     Primitive_Error("bad type");
  1276.     break;
  1277.   } /*switch*/
  1278.  
  1279.   return B;
  1280. } /*vfuncall*/
  1281.  
  1282.  
  1283. /*(DOCENTRY (USAGE "(v-sin v)"))
  1284. */
  1285. #define VSIN   Pvsin, "v-sin", 1,1,EVAL,
  1286. static Object 
  1287. Pvsin(A)  Object A;
  1288. {  Error_Tag = "v-sin";  return vfuncall(A,sin);  }
  1289.  
  1290. /*(DOCENTRY (USAGE "(v-cos v)"))
  1291. */
  1292. #define VCOS   Pvcos, "v-cos", 1,1,EVAL,
  1293. static Object 
  1294. Pvcos(A)  Object A;
  1295. {  Error_Tag = "v-cos";  return vfuncall(A,cos);  }
  1296.  
  1297. /*(DOCENTRY (USAGE "(v-sqrt v)"))
  1298. */
  1299. #define VSQRT   Pvsqrt, "v-sqrt", 1,1,EVAL,
  1300. static Object 
  1301. Pvsqrt(A)  Object A;
  1302. {  Error_Tag = "v-sqrt";  return vfuncall(A,sqrt);  }
  1303.  
  1304.  
  1305. /*(DOCENTRY (USAGE "(v-exp v)"))
  1306. */
  1307. #define VEXP   Pvexp, "v-exp", 1,1,EVAL,
  1308. static Object 
  1309. Pvexp(A)  Object A;
  1310. {  Error_Tag = "v-exp";  return vfuncall(A,exp);  }
  1311.  
  1312.  
  1313. /*(DOCENTRY (USAGE "(v-abs v)"))
  1314. */
  1315. /* generic integer/real abs function */
  1316. #define VABS   Pvabs, "v-abs", 1,1,EVAL,
  1317. Object 
  1318. Pvabs(A)
  1319.   Object A;
  1320. {
  1321.   register int i,len;
  1322.   Object B;
  1323.   Farray *a;
  1324.   GC_Node;
  1325.   Error_Tag = "v-abs";
  1326.  
  1327.   Check_Type(A,T_Farray);
  1328.  
  1329.   GC_Link(A);
  1330.   B = farray_make_like(A);
  1331.   GC_Unlink;
  1332.  
  1333.   a = FARRAY(A);
  1334.   len = a->len;
  1335.  
  1336. # define VABSTYPE(typ) \
  1337.   { \
  1338.     register typ *ia,*ib;\
  1339.     ia = (typ *)a->data;\
  1340.     ib = (typ *)FARRAY(B)->data;\
  1341.     for( i=0; i < len; i++ ) {\
  1342.       *ib++ = *ia < (typ)0 ? - *ia : *ia;\
  1343.       ia++;\
  1344.     }\
  1345.   }
  1346.  
  1347.   switch(a->type) {
  1348.  
  1349.   case T_Fixnum: VABSTYPE(int4) break;
  1350.     
  1351.   case T_Flonum: VABSTYPE(float) break;
  1352.  
  1353.   default:
  1354.     Primitive_Error("bad type");
  1355.     break;
  1356.   } /*switch*/
  1357.  
  1358.   return B;
  1359. # undef VABSTYPE
  1360. } /*vabs*/
  1361.  
  1362.  
  1363.  
  1364. /*(DOCENTRY
  1365.   (USAGE "(v-not v) -- requires and returns an boolean vector"))
  1366. */
  1367. #define VNOT   Pvnot, "v-not", 1,1,EVAL,
  1368. Object 
  1369. Pvnot(A)
  1370.   Object A;
  1371. {
  1372.   register int i,len;
  1373.   Object B;
  1374.   Farray *a;
  1375.   GC_Node;
  1376.   Error_Tag = "v-not";
  1377.  
  1378.   Check_Type(A,T_Farray);
  1379.  
  1380.   GC_Link(A);
  1381.   B = farray_make_like(A);
  1382.   GC_Unlink;
  1383.  
  1384.   a = FARRAY(A);
  1385.   len = a->len;
  1386.  
  1387.   switch(a->type) {
  1388.   case T_Vbool:  {
  1389.     register Vbool *ia,*ib;
  1390.     ia = (Vbool *)a->data;
  1391.     ib = (Vbool *)FARRAY(B)->data;
  1392.     for( i=0; i < len; i++ ) {
  1393.       *ib++ = 1 - *ia++;
  1394.     }
  1395.     break;
  1396.   }
  1397.  
  1398.   default:
  1399.     Primitive_Error("bad type");
  1400.     break;
  1401.   } /*switch*/
  1402.  
  1403.   return B;
  1404. } /*vnot*/
  1405.  
  1406.  
  1407. /*%%%%*/
  1408.  
  1409. /*(DOCENTRY
  1410.   (USAGE "(v-rnd len|array) => Returns a vector[len] of random values."
  1411.          "len may be another farray, in which case its length is used."))
  1412. */
  1413. #if ZILLAONLY
  1414. # include <rnd.h>
  1415. #else
  1416.   double rndf() { return( random() / (double)((1<<31)-1) ); }
  1417. # define rndbit() (random() & 0x1)
  1418. #endif
  1419.  
  1420. #define VRANDOM    Pvrandom, "v-rnd", 1,1,EVAL,
  1421. Object
  1422. Pvrandom(Len)
  1423.   Object Len;
  1424. {
  1425.   register int4 i,len;
  1426.   Object vd;
  1427.   Error_Tag = "v-rnd"; 
  1428.  
  1429.   if (!((TYPE(Len)==T_Fixnum)||(TYPE(Len)==T_Farray)||(TYPE(Len)==T_Bignum)))
  1430.     Primitive_Error("length must be integer or sample farray");
  1431.  
  1432.   if (TYPE(Len)==T_Farray) {
  1433.     len = FARRAY(Len)->len;
  1434.     vd = farray_make_like(Len);
  1435.  
  1436.     switch(FARRAY(vd)->type) {
  1437.     case T_Flonum: {
  1438.       FARRAYDATAPTR(iv,vd,float);
  1439.       for( i=0; i < len; i++ ) *iv++ = rndf();
  1440.     } break;
  1441.     case T_Fixnum: {
  1442.       FARRAYDATAPTR(iv,vd,int4);
  1443.       for( i=0; i < len; i++ ) *iv++ = rndbit();
  1444.     } break;
  1445.     case T_String: {
  1446.       FARRAYDATAPTR(iv,vd,char);
  1447.       for( i=0; i < len; i++ ) *iv++ = (char)rndbit();
  1448.     } break;
  1449.     }/*switch*/
  1450.   } /*len is farray*/
  1451.  
  1452.   else {
  1453.     register float *iv;
  1454.     len = Get_Integer(Len);
  1455.     vd = farray_make(T_Flonum,len);
  1456.     iv = (float *)FARRAY(vd)->data;
  1457.     for( i=0; i < len; i++ ) *iv++ = rndf();
  1458.   }
  1459.  
  1460.   return vd;
  1461. } /*vrandom*/
  1462.  
  1463.  
  1464. /*%%%%%%%%%%%%%%%% binary (dyadic) arithmetic operations %%%%%%%%%%%%%%%%*/
  1465. /*(SECTION
  1466.   "Dyadic Functions"
  1467.   "Mixed (int,float) argument types promote to float result vectors.")
  1468. */
  1469.  
  1470. /* promote mixed type expressions: int->float, string->int */
  1471. static int type_elevate(A,B)
  1472.   Object A,B;
  1473. {
  1474.   int atype,btype;
  1475.   Ztrace(("fvector type_elevate a=%x,b=%x\n",A,B));
  1476.   atype = TYPE(A); if (atype == T_Farray) atype = FARRAY(A)->type;
  1477.   btype = TYPE(B); if (btype == T_Farray) btype = FARRAY(B)->type;
  1478.  
  1479. #if 0  /* more elegant, but less control and error checking: */
  1480.   if ((atype == T_Flonum) || (btype == T_Flonum)) return T_Flonum;
  1481.   if ((atype == T_Fixnum) || (btype == T_Fixnum)) return T_Fixnum;
  1482.   Ztrace(("--fvector type_elevate\n")); fflush(stdout);
  1483.   return T_String;
  1484. #endif
  1485.  
  1486.   switch(atype) {
  1487.   case T_Flonum:
  1488.     switch(btype) {
  1489.     case T_Flonum:
  1490.     case T_Fixnum:      return T_Flonum;
  1491.     case T_String:      Primitive_Error("bad vector types: flt,string");
  1492.     default: Panic("bad array type in type_elevate");
  1493.     } break;
  1494.   case T_Fixnum:
  1495.     switch(btype) {
  1496.     case T_Flonum:      return T_Flonum;
  1497.     case T_Fixnum:      return T_Fixnum;
  1498.     case T_String:      Primitive_Error("bad vector types: fix,string");
  1499.     default: Panic("bad array type in type_elevate");
  1500.     } break;
  1501.   case T_String:
  1502.     switch(btype) {
  1503.     case T_String:      return T_String;
  1504.     default:            Primitive_Error("bad vector types: string,*");
  1505.     }
  1506.     break;
  1507.   default: Panic("bad array type in type_elevate");
  1508.   } /*switch*/
  1509.  
  1510.   Panic("type-elevate failure");
  1511. } /*type_elevate*/
  1512.  
  1513.  
  1514. /* strange function called from within ARITHOP */
  1515. static Object generic_badargs Zproto((int,Object[]));
  1516. static Object generic_badargs(ac,av) 
  1517.   int ac; Object av[];
  1518. {
  1519.   Primitive_Error("neither arg is vector");
  1520.   return Null;
  1521. }
  1522.  
  1523. #ifdef OLD /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
  1524.  
  1525. /* binary arithmetic op loop, arrayXarray */
  1526. # define ARITHBINAALOOP(op, xtype, ytype, ztype, alen) {\
  1527.     register xtype *ix = (xtype *)a->data; \
  1528.     register ytype *iy = (ytype *)b->data; \
  1529.     register ztype *iz = (ztype *)c->data; \
  1530.     register int i; \
  1531.     if ((alen%4) == 0) { /* UNROLL LOOP! */ \
  1532.       register int alen4 = alen>>2;\
  1533.       for( i=0; i < alen4; i++ ) {\
  1534.         *iz++ = op(*ix,*iy); ix++; iy++;\
  1535.         *iz++ = op(*ix,*iy); ix++; iy++;\
  1536.         *iz++ = op(*ix,*iy); ix++; iy++;\
  1537.         *iz++ = op(*ix,*iy); ix++; iy++;\
  1538.       }\
  1539.     }\
  1540.     else \
  1541.       for( i=0; i < alen; i++ ) {\
  1542.         *iz++ = op(*ix,*iy); ix++; iy++;\
  1543.       }\
  1544.   } /*arithbinaaloop*/
  1545.  
  1546.  
  1547. /* OP is operation to do when result is integer, otherwise FOP.
  1548.  * This distinction is not needed for +*-/, but it is for %,fmod()
  1549.  */
  1550. #define ARITHBINAA(AANAME,SNAME,OP,FOP) \
  1551. static Object AANAME(A,B)  \
  1552.   Object A,B; \
  1553. {\
  1554.   Object C;\
  1555.   Farray *a,*b,*c;\
  1556.   register int len;\
  1557.   GC_Node2;\
  1558.   Error_Tag = SNAME;\
  1559.                 \
  1560.   len = v_conform(A,B);\
  1561.   a = FARRAY(A); b = FARRAY(B);\
  1562.   Ztrace(("%s a[%d]type=%d,  b[%d]type=%d\n",\
  1563.           SNAME,a->len,a->type,b->len,b->type));\
  1564.                 \
  1565.   GC_Link2(A,B);\
  1566.   C = farray_make(type_elevate(A,B),b->len); \
  1567.   GC_Unlink;\
  1568.   a = FARRAY(A); b = FARRAY(B);  /*reassign after gc*/\
  1569.   c = FARRAY(C);\
  1570.                 \
  1571.   /* four cases: int*int, flt*int, int*flt, flt*flt */ \
  1572.                 \
  1573.   switch(a->type) {\
  1574.   case T_Fixnum:\
  1575.     if (b->type == T_Fixnum) {\
  1576.       ARITHBINAALOOP(OP, int4,int4,int4, len)\
  1577.       break;\
  1578.     }\
  1579.     else if (b->type == T_Flonum) {\
  1580.       ARITHBINAALOOP(FOP, int4,float, float, len)\
  1581.       break;\
  1582.     }\
  1583.     else goto err;\
  1584.     break;\
  1585.                 \
  1586.   case T_Flonum: \
  1587.     if (b->type == T_Flonum) {\
  1588.       ARITHBINAALOOP(FOP, float,float, float, len)\
  1589.       break;\
  1590.     }\
  1591.     else if (b->type == T_Fixnum) {\
  1592.       ARITHBINAALOOP(FOP, float,int4,float, len)\
  1593.       break;\
  1594.     }\
  1595.     else goto err;\
  1596.     break;\
  1597.                 \
  1598.   default: goto err;\
  1599.   } /*switch*/\
  1600.                 \
  1601.   return C;\
  1602.   err:;\
  1603.   Primitive_Error("bad vector types");\
  1604. } /*ARITHBINAA*/
  1605.  
  1606. #endif /*OLD%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
  1607.  
  1608.  
  1609. /* binary arithmetic op loop, arrayXarray */
  1610. # define ARITHOPLOOP(op,ztype,alen, xinit,x,xinc, yinit,y,yinc) {\
  1611.     xinit \
  1612.     yinit \
  1613.     register ztype *iz = (ztype *)c->data; \
  1614.     register int i; \
  1615.     if ((alen%4) == 0) { /* UNROLL LOOP! */ \
  1616.       register int alen4 = alen>>2;\
  1617.       for( i=0; i < alen4; i++ ) {\
  1618.         *iz++ = op(x,y); xinc yinc \
  1619.         *iz++ = op(x,y); xinc yinc \
  1620.         *iz++ = op(x,y); xinc yinc \
  1621.         *iz++ = op(x,y); xinc yinc \
  1622.       }\
  1623.     }\
  1624.     else \
  1625.       for( i=0; i < alen; i++ ) {\
  1626.         *iz++ = op(x,y); xinc yinc \
  1627.       }\
  1628.   } /*arithoploop*/
  1629.  
  1630. #define A_ARRAYINT register int4 *ix = FARRAYDATA(A,int4);
  1631. #define B_ARRAYINT register int4 *iy = FARRAYDATA(B,int4);
  1632. #define A_ARRAYFLT register float *ix = FARRAYDATA(A,float);
  1633. #define B_ARRAYFLT register float *iy = FARRAYDATA(B,float);
  1634. #define A_AREF *ix
  1635. #define A_AINC ix++;
  1636. #define B_AREF *iy
  1637. #define B_AINC iy++;
  1638.  
  1639. #define A_SCALARINT register int4 x = Get_Integer(A);
  1640. #define B_SCALARINT register int4 y = Get_Integer(B);
  1641.  
  1642. #define A_SCALARFLT register float x = FLONUM(A)->val;
  1643. #define B_SCALARFLT register float y = FLONUM(B)->val;
  1644.  
  1645. #define ARITHCASE(OP,atype,btype,rtype) \
  1646.       if (aisarray && bisarray)\
  1647.         ARITHOPLOOP(OP,rtype,len,\
  1648.                     A_ARRAY##atype,A_AREF,A_AINC,\
  1649.                     B_ARRAY##btype,B_AREF,B_AINC)\
  1650.       else if (aisarray) \
  1651.         ARITHOPLOOP(OP,rtype,len,\
  1652.                     A_ARRAY##atype,A_AREF,A_AINC,\
  1653.                     B_SCALAR##btype,y,)\
  1654.       else if (bisarray) \
  1655.         ARITHOPLOOP(OP,rtype,len,\
  1656.                     A_SCALAR##atype,x,,\
  1657.                     B_ARRAY##btype,B_AREF,B_AINC)
  1658.  
  1659.  
  1660. #define ARITHOP(AANAME,SNAME,OP,ELSESTATEMENT) \
  1661. static Object AANAME(A,B)  \
  1662.   Object A,B; \
  1663. {\
  1664.   Object C;\
  1665.   Farray *a,*b,*c;\
  1666.   int atype,btype;\
  1667.   bool aisarray,bisarray;\
  1668.   register int len;\
  1669.   GC_Node2;\
  1670.   Error_Tag = SNAME;\
  1671.   Ztrace(("fvector arithop %s\n", # SNAME));\
  1672.                 \
  1673.   atype = TYPEORFARRAYTYPE(A);\
  1674.   btype = TYPEORFARRAYTYPE(B);\
  1675.   aisarray = TYPE(A)==T_Farray;\
  1676.   bisarray = TYPE(B)==T_Farray;\
  1677.   if (aisarray && bisarray) \
  1678.     len = v_conform(A,B);\
  1679.   else if (aisarray) len = FARRAY(A)->len; \
  1680.   else if (bisarray) len = FARRAY(B)->len; \
  1681.   if (!aisarray && !bisarray) {\
  1682.     Object av[2]; \
  1683.     av[0] = A; av[1] = B; \
  1684.     /* expands to either e.g. P_Generic_Minus or generic_badargs */\
  1685.     return ELSESTATEMENT (2,av);\
  1686.   } \
  1687.   GC_Link2(A,B);\
  1688.   C = v_elevated_array(A,B);\
  1689.   GC_Unlink;\
  1690.   c = FARRAY(C);\
  1691.                 \
  1692.   /* cases: int*int, flt*int, int*flt, flt*flt */ \
  1693.                 \
  1694.   switch(atype) {\
  1695.   case T_Fixnum:\
  1696.     if (btype == T_Fixnum) {\
  1697.       ARITHCASE(OP,INT,INT,int4) \
  1698.       break;\
  1699.     }\
  1700.     else if (btype == T_Flonum) {\
  1701.       ARITHCASE(OP,INT,FLT,float) \
  1702.       break;\
  1703.     }\
  1704.     else goto err;\
  1705.     break;\
  1706.                 \
  1707.   case T_Flonum: \
  1708.     if (btype == T_Fixnum) {\
  1709.       ARITHCASE(OP,FLT,INT,float) \
  1710.       break;\
  1711.     }\
  1712.     else if (btype == T_Flonum) {\
  1713.       ARITHCASE(OP,FLT,FLT,float) \
  1714.       break;\
  1715.     }\
  1716.     else goto err;\
  1717.     break;\
  1718.                 \
  1719.   default: goto err;\
  1720.   } /*switch*/\
  1721.                 \
  1722.   return C;\
  1723.   err:;\
  1724.   Primitive_Error("bad vector types");\
  1725.   return C; /*for lint*/\
  1726. } /*ARITHOP*/
  1727.  
  1728. typedef Object (elkgenericop_t) Zproto((int,Object []));
  1729. elkgenericop_t P_Generic_Plus,P_Generic_Minus,
  1730.   P_Generic_Multiply,P_Generic_Divide;
  1731.  
  1732.  
  1733. /*(DOCENTRY
  1734.   (USAGE "(v-+ a b) --  elementwise addition"))
  1735. */
  1736. #define VADD    P_vadd, "v-+", 2,2,EVAL,
  1737. #define plus(a,b) (a)+(b)
  1738. ARITHOP(P_vadd,"v-+",plus,P_Generic_Plus)
  1739. #undef plus
  1740.  
  1741. /*(DOCENTRY
  1742.   (USAGE "(v-- a b) --  elementwise subtraction"))
  1743. */
  1744. #define VSUB    P_vsub, "v--", 2,2,EVAL,
  1745. #define sub(a,b) (a)-(b)
  1746. ARITHOP(P_vsub,"v--",sub,P_Generic_Minus)
  1747. #undef sub
  1748.  
  1749.  
  1750. /*(DOCENTRY
  1751.   (USAGE "(v-* a b) --  elementwise multiplication"))
  1752. */
  1753. #define VMUL    P_vmul, "v-*", 2,2,EVAL,
  1754. #define mul(a,b) (a)*(b)
  1755. ARITHOP(P_vmul,"v-*",mul,P_Generic_Multiply)
  1756. #undef mul
  1757.  
  1758. /*(DOCENTRY
  1759.   (USAGE "(v-/ a b) --  elementwise division."
  1760.          "There is NO divide-by-zero check."))
  1761. */
  1762. #define VDIV    P_vdiv, "v-/", 2,2,EVAL,
  1763. #define div(a,b) (a)/(b)
  1764. ARITHOP(P_vdiv,"v-/",div,P_Generic_Divide)
  1765. #undef div
  1766.  
  1767. /*(DOCENTRY
  1768.   (USAGE "(v-min a b) -- elementwise min"))
  1769. */
  1770. #define VMIN    P_vmin, "v-min", 2,2,EVAL,
  1771. #define min(a,b) ((a)<(b) ? (a) : (b))
  1772. ARITHOP(P_vmin,"v-min",min,generic_badargs)
  1773. #undef min
  1774.  
  1775. /*(DOCENTRY
  1776.   (USAGE "(v-max a b) -- elementwise max"))
  1777. */
  1778. #define VMAX    P_vmax, "v-max", 2,2,EVAL,
  1779. #define max(a,b) ((a)>(b) ? (a) : (b))
  1780. ARITHOP(P_vmax,"v-max",max,generic_badargs)
  1781. #undef max
  1782.  
  1783. /*(DOCENTRY
  1784.   (USAGE "(v-mod a b) -- elementwise mod, implemented as fmod()"))
  1785. */
  1786. #define VMOD    P_vmod, "v-mod", 2,2,EVAL,
  1787. #define _fmod(a,b) fmod((double)a,(double)b)
  1788. ARITHOP(P_vmod,"v-mod",_fmod,generic_badargs)
  1789. #undef _fmod
  1790.  
  1791.  
  1792. #if ZILLAONLY
  1793. #define VTEST    junktest, "v-test", 1,1,EVAL,
  1794. static Object junktest(Object B)
  1795. {
  1796.   Error_Tag = "v-test";
  1797.  
  1798.   if (TYPE(B)==T_Farray) {
  1799.     fprintf(stderr,"B is farray\n");
  1800.     fprintf(stderr,"B->len %d",FARRAY(B)->len);
  1801.   }
  1802.   else if (TYPE(B)==T_Flonum) {
  1803.     float f;
  1804.     int adr;
  1805.     adr = (int)(&FLONUM(B)->val);
  1806.     fprintf(stderr,"B is float\n");
  1807.     fprintf(stderr,"B at %d\n",B);
  1808.     fprintf(stderr,"val ptr @%d, 4=%d, 8=%d\n",
  1809.             adr,Zalignedon((char *)adr,4), Zalignedon((char *)adr,8));
  1810.     f = (float)(FLONUM(B)->val);
  1811.     fprintf(stderr,"B->flt %f",f);
  1812.   }
  1813.   else if (TYPE(B)==T_Fixnum) {
  1814.     fprintf(stderr,"B is fix\n");
  1815.     fprintf(stderr,"B at %d\n",B);
  1816.     fprintf(stderr,"B->fix :%d:",Get_Integer(B));
  1817.   }
  1818.  
  1819.   else  Primitive_Error("unknown type for B");
  1820. } /*junk-test*/
  1821. #else
  1822. # define VTEST /*nothing*/
  1823. #endif /*ZILLAONLY*/
  1824.  
  1825. #undef ARITHCASE
  1826. #undef ARITHOPLOOP
  1827. #undef ARITHOP
  1828.  
  1829. #undef A_ARRAYINT
  1830. #undef B_ARRAYINT
  1831. #undef A_ARRAYFLT
  1832. #undef B_ARRAYFLT
  1833. #undef A_AREF
  1834. #undef A_AINC
  1835. #undef B_AREF
  1836. #undef B_AINC
  1837.  
  1838. #undef A_SCALARINT
  1839. #undef B_SCALARINT
  1840.  
  1841. #undef A_SCALARFLT
  1842. #undef B_SCALARFLT
  1843.  
  1844.  
  1845. /*(DOCENTRY
  1846.   (USAGE "(v-imod a b) -- elementwise mod, implemented as %"))
  1847. */
  1848. #define VIMOD    P_vimod, "v-imod", 2,2,EVAL,
  1849.  
  1850. static Object P_vimod(A,B)
  1851.   Object A,B;
  1852. {
  1853.   Object C;
  1854.   Farray *a,*b,*c;
  1855.   register int4 *ia,*ib,*ic;
  1856.   register int i,len;
  1857.   GC_Node2;
  1858.   Error_Tag = "v-imod";
  1859.  
  1860.   len = v_conform(A,B);
  1861.   a = FARRAY(A); b = FARRAY(B);
  1862.   Ztrace(("v-binfunc a[%d]type=%d, b[%d]type=%d\n",
  1863.           a->len,a->type,b->len,b->type));
  1864.  
  1865.   if ((a->type != T_Fixnum) || (b->type != T_Fixnum))
  1866.     Primitive_Error("need integer farrays");
  1867.  
  1868.   GC_Link2(A,B);
  1869.   C = farray_make_like(A);
  1870.   GC_Unlink;
  1871.  
  1872.   a = FARRAY(A); b = FARRAY(B);
  1873.   c = FARRAY(C);
  1874.  
  1875.   ia = (int4 *)a->data; ib = (int4 *)b->data; ic = (int4 *)c->data;
  1876.   for( i=0; i < len; i++ )  {
  1877.     *ic++ = *ia % *ib;
  1878.     ia++; ib++;
  1879.   }
  1880.  
  1881.   return C;
  1882. } /*vimod*/
  1883.  
  1884.  
  1885.  
  1886. /*%%%%%%%%%%%%%%%% boolean ops %%%%%%%%%%%%%%%%*/
  1887. /*(SECTION
  1888.   "Boolean Functions"
  1889.   "The result type is an integer vector.")
  1890. */
  1891.  
  1892. # define BLOOP(op, xtype, ytype, len) { \
  1893.   register xtype *ix = (xtype *)a->data; \
  1894.   register ytype *iy = (ytype *)b->data; \
  1895.   register Vbool *iz = (Vbool *)c->data; \
  1896.   register int i; \
  1897.   for( i=0; i < len; i++ )  {\
  1898.     *iz++ = op(*ix,*iy); ix++; iy++;\
  1899.   }}
  1900.  
  1901. #define VBOP(NAME,SNAME,OP) \
  1902. Object NAME(A,B)  \
  1903.   Object A,B; \
  1904. {\
  1905.   Object C;\
  1906.   Farray *a,*b,*c;\
  1907.   register int len;\
  1908.   GC_Node2;\
  1909.   Error_Tag = SNAME;\
  1910.                 \
  1911.   len = v_conform(A,B);\
  1912.   a = FARRAY(A); b = FARRAY(B);\
  1913.   Ztrace(("%s a[%d]type=%d,  b[%d]type=%d\n",\
  1914.           SNAME,a->len,a->type,b->len,b->type));\
  1915.                 \
  1916.   GC_Link2(A,B);\
  1917.   C = farray_make(T_Vbool,b->len);\
  1918.   GC_Unlink;\
  1919.   farray_copyshape(A,C);\
  1920.     \
  1921.   a = FARRAY(A); b = FARRAY(B); /*reassign after gc*/\
  1922.   c = FARRAY(C);\
  1923.                 \
  1924.   /* four cases: int*int, flt*int, int*flt, flt*flt */ \
  1925.                 \
  1926.   switch(a->type) {\
  1927.   case T_Fixnum:\
  1928.     if (b->type == T_Fixnum) {\
  1929.       BLOOP(OP, int4, int4, len)\
  1930.       break;\
  1931.     }\
  1932.     else if (b->type == T_Flonum) {\
  1933.       BLOOP(OP, int4, float, len)\
  1934.       break;\
  1935.     }\
  1936.     else goto err;\
  1937.     break;\
  1938.                 \
  1939.   case T_Flonum: \
  1940.     if (b->type == T_Flonum) {\
  1941.       BLOOP(OP, float, float, len)\
  1942.       break;\
  1943.     }\
  1944.     else if (b->type == T_Fixnum) {\
  1945.       BLOOP(OP, float, int4, len)\
  1946.       break;\
  1947.     }\
  1948.     else goto err;\
  1949.     break;\
  1950.                 \
  1951.   default: goto err;\
  1952.   } /*switch*/\
  1953.                 \
  1954.   return C;\
  1955.                 \
  1956.   err:;\
  1957.   Primitive_Error("bad vector types");\
  1958.   return C; /*for lint*/\
  1959. } /*vbop*/
  1960.  
  1961.  
  1962. /*(DOCENTRY (USAGE "(v-eq a b)"))
  1963. */
  1964. #define VEQ    Pveq, "v-eq", 2,2,EVAL,
  1965. #define _eq(a,b) (a)==(b)
  1966. VBOP(Pveq,"v-eq",_eq)
  1967. #undef _eq
  1968.  
  1969. /*(DOCENTRY (USAGE "(v-ne a b)"))
  1970. */
  1971. #define VNE    Pvne, "v-ne", 2,2,EVAL,
  1972. #define _ne(a,b) (a)!=(b)
  1973. VBOP(Pvne,"v-ne",_ne)
  1974. #undef _ne
  1975.  
  1976. /*(DOCENTRY (USAGE "(v-lt a b)"))
  1977. */
  1978. #define VLT    Pvlt, "v-lt", 2,2,EVAL,
  1979. #define _lt(a,b) (a)<(b)
  1980. VBOP(Pvlt,"v-lt",_lt)
  1981. #undef _lt
  1982.  
  1983. /*(DOCENTRY (USAGE "(v-le a b)"))
  1984. */
  1985. #define VLE    Pvle, "v-le", 2,2,EVAL,
  1986. #define _le(a,b) (a)<=(b)
  1987. VBOP(Pvle,"v-le",_le)
  1988. #undef _le
  1989.  
  1990. /*(DOCENTRY (USAGE "(v-gt a b)"))
  1991. */
  1992. #define VGT    Pvgt, "v-gt", 2,2,EVAL,
  1993. #define _gt(a,b) (a)>(b)
  1994. VBOP(Pvgt,"v-gt",_gt)
  1995. #undef _gt
  1996.  
  1997. /*(DOCENTRY (USAGE "(v-ge a b)"))
  1998. */
  1999. #define VGE    Pvge, "v-ge", 2,2,EVAL,
  2000. #define _ge(a,b) (a)>=(b)
  2001. VBOP(Pvge,"v-ge",_ge)
  2002. #undef _ge
  2003.  
  2004. /*(DOCENTRY (USAGE "(v-and a b)"))
  2005. */
  2006. #define VAND    Pvand, "v-and", 2,2,EVAL,
  2007. #define _and(a,b) (a)&&(b)
  2008. VBOP(Pvand,"v-and",_and)
  2009. #undef _and
  2010.  
  2011. /*(DOCENTRY (USAGE "(v-or a b)"))
  2012. */
  2013. #define VOR    Pvor, "v-or", 2,2,EVAL,
  2014. #define _or(a,b) (a)||(b)
  2015. VBOP(Pvor,"v-or",_or)
  2016. #undef _or
  2017.  
  2018. /*(DOCENTRY
  2019.   (USAGE "(v-xor a b) -- ((a||b)&&(!(a&&b)))"))
  2020. */
  2021. #define VXOR    Pvxor, "v-xor", 2,2,EVAL,
  2022. #define _xor(a,b) (((a)||(b)) && !((a)&&(b)))
  2023. VBOP(Pvxor,"v-xor",_xor)
  2024. #undef _xor
  2025.  
  2026. #undef BLOOP
  2027. #undef VBOP
  2028.  
  2029.  
  2030. /*%%%%%%%%%%%%%%%% float dyadic functions %%%%%%%%%%%%%%%%*/
  2031. /*(SECTION
  2032.   "Float Dyadic Functions"
  2033.   "These functions always return float vectors,"
  2034.   "rather than vectors of the highest argument type.")
  2035. */
  2036.  
  2037. static Object vbinfunc
  2038.     Zproto((Object,Object,double(*)(double,double),char *));
  2039. static Object vbinfunc(A,B,f,name)
  2040.   Object A,B;
  2041.   double (*f) Zproto((double,double));
  2042.   char *name;
  2043. {
  2044.   Object C;
  2045.   Farray *a,*b,*c;
  2046.   register float *ia,*ib,*ic;
  2047.   register int i,len;
  2048.   GC_Node2;
  2049.   Error_Tag = name;
  2050.  
  2051.   len = v_conform(A,B);
  2052.   /* may not be commutative, do not reorder! */
  2053.   a = FARRAY(A); b = FARRAY(B);
  2054.   Ztrace(("v-binfunc a[%d]type=%d, b[%d]type=%d\n",
  2055.           a->len,a->type,b->len,b->type));
  2056.  
  2057.   if ((a->type != T_Flonum) || (b->type != T_Flonum))
  2058.     Primitive_Error("need float farrays");
  2059.  
  2060.   GC_Link2(A,B);
  2061.   C = farray_make(T_Flonum,b->len);
  2062.   GC_Unlink;
  2063.   farray_copyshape(A,C);
  2064.  
  2065.   a = FARRAY(A); b = FARRAY(B);
  2066.   c = FARRAY(C);
  2067.  
  2068.   ia = (float *)a->data; ib = (float *)b->data; ic = (float *)c->data;
  2069.   for( i=0; i < len; i++ )  {
  2070.     *ic++ = (*f)(*ia,*ib);
  2071.     ia++; ib++;
  2072.   }
  2073.  
  2074.   return C;
  2075. } /*vbinfunc*/
  2076.  
  2077.  
  2078. /*(DOCENTRY
  2079.   (USAGE "(v-pow x p) -- elementwise power function"))
  2080. */
  2081. #define VPOW    Pvpow, "v-pow", 2,2,EVAL,
  2082. Object Pvpow(a,b) Object a,b; { return vbinfunc(a,b,pow,"v-pow"); }
  2083.  
  2084. /*(DOCENTRY
  2085.   (USAGE "(v-atan2 y x) -- elementwise atan2(y,x)."))
  2086. */
  2087. #define VATAN2    Pvatan2, "v-atan2", 2,2,EVAL,
  2088. Object Pvatan2(a,b) Object a,b; { return vbinfunc(a,b,atan2,"v-atan2"); }
  2089.  
  2090.  
  2091. /*%%%%%%%%%%%%%%%% REDUCTIONS %%%%%%%%%%%%%%%%*/
  2092. /*(SECTION
  2093.   "Reductions"
  2094.   "scalar = (op-reduce arr)"
  2095.   "---------------"
  2096.   "scalar r = arr[0];"
  2097.   "for( i=1; i < len(arr); i++ ) r = op(r,arr[i]);")
  2098. */
  2099.  
  2100. # define REDUCELOOP(a,type,len, op) {\
  2101.     register type *arr = (type *)a->data; \
  2102.     register int i; \
  2103.     z = *arr++; /*z=arr[0]. (z is free variable) */\
  2104.     for( i=1; i < len; i++ ) {\
  2105.       z = op(z,*arr); arr++;\
  2106.     }\
  2107.   }
  2108.  
  2109. /* gc protect not necessary here */
  2110. #define REDUCE(NAME,SNAME,OP) \
  2111. static Object NAME(A)  \
  2112.   Object A; \
  2113. {\
  2114.   Object B;\
  2115.   Farray *a;\
  2116.   register int len;\
  2117.   Error_Tag = SNAME;\
  2118.   \
  2119.   Check_Type(A,T_Farray);\
  2120.   a = FARRAY(A); len = a->len;\
  2121.   Ztrace(("%s a[%d]type=%d\n",SNAME,a->len,a->type));\
  2122.                 \
  2123.   switch(a->type) {\
  2124.   case T_Fixnum: {\
  2125.     int4 z;\
  2126.     REDUCELOOP(a,int4,len,OP)\
  2127.     B = Make_Integer(z);\
  2128.     break;\
  2129.   }\
  2130.   case T_Flonum: {\
  2131.     float z;\
  2132.     REDUCELOOP(a,float,len,OP)\
  2133.     B = Make_Reduced_Flonum(z);\
  2134.     break;\
  2135.   }\
  2136.   case T_String: {\
  2137.     Zbyte z;\
  2138.     REDUCELOOP(a,Zbyte,len,OP)\
  2139.     B = Make_Integer((int4)z);\
  2140.     break;\
  2141.   }\
  2142.   default: goto err;\
  2143.   } /*switch*/\
  2144.                 \
  2145.   return B;\
  2146.                 \
  2147.   err:;\
  2148.   Primitive_Error("bad vector types");\
  2149.   return B; /*for lint*/\
  2150. } /*REDUCE*/
  2151.  
  2152.  
  2153. /*(DOCENTRY (USAGE "(+/ v) or v-+reduce"))
  2154. */
  2155. #define VPLUSREDUCE    P_vplusreduce, "v-+reduce", 1,1,EVAL,
  2156. #define VPLUSREDUCEb   P_vplusreduce, "v-+/", 1,1,EVAL,
  2157. #define VPLUSREDUCEc   P_vplusreduce, "+/", 1,1,EVAL,
  2158. #define _plus(a,b) ((a)+(b))
  2159. REDUCE(P_vplusreduce,"vplusreduce",_plus)
  2160. #undef _plus
  2161.  
  2162.  
  2163. /*(DOCENTRY (USAGE "(-/ v) or v--reduce"))
  2164. */
  2165. #define VMINUSREDUCE    P_vminusreduce, "v--reduce", 1,1,EVAL,
  2166. #define VMINUSREDUCEb   P_vminusreduce, "v--/", 1,1,EVAL,
  2167. #define VMINUSREDUCEc   P_vminusreduce, "-/", 1,1,EVAL,
  2168. #define _minus(a,b) ((a)-(b))
  2169. REDUCE(P_vminusreduce,"vminusreduce",_minus)
  2170. #undef _minus
  2171.  
  2172. /*(DOCENTRY (USAGE "(* / v) or v-*reduce"))
  2173. */
  2174. #define VMULREDUCE    P_vmulreduce, "v-*reduce", 1,1,EVAL,
  2175. #define VMULREDUCEb   P_vmulreduce, "v-*/", 1,1,EVAL,
  2176. #define VMULREDUCEc   P_vmulreduce, "*/", 1,1,EVAL,
  2177. #define _mul(a,b) ((a)*(b))
  2178. REDUCE(P_vmulreduce,"mulreduce",_mul)
  2179. #undef _mul
  2180.  
  2181. /*(DOCENTRY (USAGE "(min/ v) or v-minreduce"))
  2182. */
  2183. #define VMINREDUCE    P_vminreduce, "v-minreduce", 1,1,EVAL,
  2184. #define VMINREDUCEb   P_vminreduce, "v-min/", 1,1,EVAL,
  2185. #define VMINREDUCEc   P_vminreduce, "min/", 1,1,EVAL,
  2186. #define _min(a,b) ((a)<=(b) ? (a) : (b))
  2187. REDUCE(P_vminreduce,"minreduce",_min)
  2188. #undef _min
  2189.  
  2190. /*(DOCENTRY (USAGE "(max/ v) or v-maxreduce"))
  2191. */
  2192. #define VMAXREDUCE    P_vmaxreduce, "v-maxreduce", 1,1,EVAL,
  2193. #define VMAXREDUCEb   P_vmaxreduce, "v-max/", 1,1,EVAL,
  2194. #define VMAXREDUCEc   P_vmaxreduce, "max/", 1,1,EVAL,
  2195. #define _max(a,b) ((a)>=(b) ? (a) : (b))
  2196. REDUCE(P_vmaxreduce,"maxreduce",_max)
  2197. #undef _max
  2198.  
  2199. #undef REDUCELOOP
  2200. #undef REDUCE
  2201.  
  2202.  
  2203. /*%%%%%%%%%%%%%%%% scan operators %%%%%%%%%%%%%%%%*/
  2204. /*(SECTION
  2205.   "Scans"
  2206.   "outvector[0] = z = vector[0];"
  2207.   "for( i=1; i < len(arr); i++ ) {"
  2208.   "  z = op(z,vector[i]);"
  2209.   "  outvector[i] = z;"
  2210.   "}")
  2211. */
  2212.  
  2213. /* UNLIKE BLELLOCH:
  2214.  * blelloch +-scan (forexample) is define as resulting in [0,...],
  2215.  * which means that the last element of the input vector is IGNORED,
  2216.  * and a non-informative zero is prepended.
  2217.  * Instead, shift so all input numbers are reflected.
  2218.  * we can make a zero anytime.
  2219.  *
  2220.  * a drawback is that the index operation (make vector of index vals)
  2221.  * cannot then be defined as (+-scan (distribute 1 len))
  2222.  */
  2223.  
  2224.  
  2225. # define SCANLOOP(arr,type,len,rslt, op) { \
  2226.     register type *ap = (type *)arr->data; \
  2227.     register type *r = (type *)rslt->data; \
  2228.     register int i; \
  2229.     register type z;\
  2230.     *r++ = z = *ap++;     /* z = ap[0] */\
  2231.     for( i=1; i < len; i++ )  {\
  2232.       *r++ = z = op(z,*ap); ap++;\
  2233.     }\
  2234.   }
  2235.  
  2236. #define SCAN(NAME,SNAME,OP) \
  2237. static Object NAME(A)  \
  2238.   Object A; \
  2239. {\
  2240.   Object B;\
  2241.   Farray *a,*b;\
  2242.   register int len;\
  2243.   GC_Node;\
  2244.   Error_Tag = SNAME;\
  2245.                 \
  2246.   Ztrace(("%s a[%d]type=%d\n",SNAME,FARRAY(A)->len,FARRAY(A)->type));\
  2247.   Check_Type(A,T_Farray); \
  2248.   len = FARRAY(A)->len;\
  2249.   GC_Link(A);\
  2250.   B = farray_make_like(A);\
  2251.   GC_Unlink;\
  2252.   a = FARRAY(A); b = FARRAY(B);\
  2253.                 \
  2254.   switch(a->type) {\
  2255.   case T_Fixnum:\
  2256.     SCANLOOP(a,int4,len,b, OP)\
  2257.     break;\
  2258.   case T_Flonum:\
  2259.     SCANLOOP(a,float,len,b, OP)\
  2260.     break;\
  2261.   case T_String:\
  2262.     SCANLOOP(a,Zbyte,len,b, OP)\
  2263.     break;\
  2264.   default: goto err;\
  2265.   } /*switch*/\
  2266.                 \
  2267.   return B;\
  2268.                 \
  2269.   err:;\
  2270.   Primitive_Error("bad vector types");\
  2271.   return B; /*for lint*/\
  2272. } /*SCAN*/
  2273.  
  2274.  
  2275. /*(DOCENTRY (USAGE "(\\+ v) or v-+scan"))
  2276. */
  2277. #define VPLUSSCAN    P_vplusscan, "v-+scan", 1,1,EVAL,
  2278. #define VPLUSSCANc   P_vplusscan, "\\+", 1,1,EVAL,
  2279. #define _plus(a,b) ((a)+(b))
  2280. SCAN(P_vplusscan,"vplusscan",_plus)
  2281. #undef _plus
  2282.  
  2283.  
  2284. /*(DOCENTRY (USAGE "(\\- v) or v--scan"))
  2285. */
  2286. #define VMINUSSCAN    P_vminusscan, "v--scan", 1,1,EVAL,
  2287. #define VMINUSSCANc   P_vminusscan, "\\-", 1,1,EVAL,
  2288. #define _minus(a,b) ((a)-(b))
  2289. SCAN(P_vminusscan,"vminusscan",_minus)
  2290. #undef _minus
  2291.  
  2292. /*(DOCENTRY (USAGE "(\\* v) or v-*scan"))
  2293. */
  2294. #define VMULSCAN    P_vmulscan, "v-*scan", 1,1,EVAL,
  2295. #define VMULSCANc   P_vmulscan, "\\*", 1,1,EVAL,
  2296. #define _mul(a,b) ((a)*(b))
  2297. SCAN(P_vmulscan,"mulscan",_mul)
  2298. #undef _mul
  2299.  
  2300. /*(DOCENTRY (USAGE "(\\min v) or v-minscan"))
  2301. */
  2302. #define VMINSCAN    P_vminscan, "v-minscan", 1,1,EVAL,
  2303. #define VMINSCANc   P_vminscan, "\\min", 1,1,EVAL,
  2304. #define _min(a,b) ((a)<=(b) ? (a) : (b))
  2305. SCAN(P_vminscan,"minscan",_min)
  2306. #undef _min
  2307.  
  2308. /*(DOCENTRY (USAGE "(\\max v) or v-maxscan"))
  2309. */
  2310. #define VMAXSCAN    P_vmaxscan, "v-maxscan", 1,1,EVAL,
  2311. #define VMAXSCANc   P_vmaxscan, "\\max", 1,1,EVAL,
  2312. #define _max(a,b) ((a)>=(b) ? (a) : (b))
  2313. SCAN(P_vmaxscan,"maxscan",_max)
  2314. #undef _max
  2315.  
  2316. #undef SCANLOOP
  2317. #undef SCAN
  2318.  
  2319. /*%%%%%%%%%%%%%%%% distribute, index %%%%%%%%%%%%%%%%*/
  2320. /*(SECTION "Miscellaneous")
  2321. */
  2322.  
  2323. /*(MANENTRY
  2324.   "(v-distribute value len)"
  2325.   "Generate a vector[len] initialized to value.")
  2326. */
  2327. #define VDISTRIBUTE    Pvdistribute, "v-distribute", 2,2,EVAL,
  2328. Object
  2329. Pvdistribute(Value,Len)
  2330.   Object Value;
  2331.   Object Len;
  2332. {
  2333.   register int4 i,len;
  2334.   Object vd;
  2335.   Error_Tag = "vdistribute";
  2336.   /* no need to gc_link as long as both len,value are retrieved
  2337.      before array is created */
  2338.  
  2339.   if (!((TYPE(Len)==T_Fixnum)||(TYPE(Len)==T_Farray)||(TYPE(Len)==T_Bignum)))
  2340.     Primitive_Error("length must be integer or sample farray");
  2341.  
  2342.   if (TYPE(Len)==T_Farray)
  2343.     len = FARRAY(Len)->len;
  2344.   else
  2345.     len = Get_Integer(Len);
  2346.  
  2347.   switch(TYPE(Value)) {
  2348.  
  2349.   case T_Fixnum: {
  2350.     int4 val = Get_Integer(Value);
  2351.     int4 *iv;
  2352.     vd = farray_make(T_Fixnum,len);
  2353.     iv = (int4 *)(FARRAY(vd)->data);
  2354.     for( i=0; i < len; i++ ) *iv++ = val;
  2355.     break;
  2356.   }
  2357.  
  2358.   case T_Flonum: {
  2359.     float val = (float)FLONUM(Value)->val;
  2360.     float *iv;
  2361.     vd = farray_make(T_Flonum,len);
  2362.     iv = (float *)(FARRAY(vd)->data);
  2363.     for( i=0; i < len; i++ ) *iv++ = val;
  2364.     break;
  2365.   }
  2366.  
  2367.   default: Primitive_Error("bad type");
  2368.   } /*switch*/
  2369.  
  2370.   return( vd );
  2371. } /*vdistribute*/
  2372.  
  2373.  
  2374. /*(MANENTRY
  2375.   "(v-emote v)"
  2376.   "Squeeze emotional value from a vector.")
  2377. */
  2378.  
  2379.  
  2380. /*(MANENTRY
  2381.   "(v-index len)"
  2382.   "Generate a vector[length] of index values."
  2383.   "(v-index 3) => (% 0 1 2)")
  2384. */
  2385. /* see blelloch p.65 */
  2386. #define VINDEX    Pvindex, "v-index", 1,1,EVAL,
  2387. Object
  2388. Pvindex(Len)
  2389.   Object Len;
  2390. {
  2391.   register int4 i,len;
  2392.   Object vd;
  2393.   register int4 *iv;
  2394.   Error_Tag = "v-index";
  2395.  
  2396.   len = Get_Integer(Len);
  2397.   Ztrace(("v-index [%d]\n",len));
  2398.  
  2399.   vd = farray_make(T_Fixnum,len);
  2400.   iv = (int4 *)FARRAY(vd)->data;
  2401.  
  2402.   for( i=0; i < len; i++ ) *iv++ = i;
  2403.  
  2404.   return vd;
  2405. } /*vindex*/
  2406.  
  2407.  
  2408. /* LOOPFUNC(name,inita,initb,statement)
  2409.    Generates a function B = name(A).
  2410.    Statement is '*bp++ = *ap++' to simply copy.
  2411.    Inita,initb are inserted after ap=A->data, bp=B->data.
  2412.    To reverse A,
  2413.      LOOPFUNC(Pvreverse,,bp += (len-1);,*bp-- = *ap++;)
  2414.    The statements inita,initb,statement should include their
  2415.    trailing semicolons.
  2416.  */
  2417.  
  2418. /* helper to LOOPFUNC, also used separately */
  2419. # define VECLOOP(atype,btype,inita,initb,statement) {\
  2420.     atype *ap; btype *bp;\
  2421.     ap = (atype *)a->data;\
  2422.     bp = (btype *)b->data;\
  2423.     inita initb\
  2424.     for( i=0; i < len; i++ ) {\
  2425.       statement\
  2426.     }\
  2427.   }
  2428.  
  2429.  
  2430. #define LOOPFUNC(name,inita,initb,statement) \
  2431. static Object name(A) \
  2432.   Object A; \
  2433. { \
  2434.   Object B;\
  2435.   Farray *a,*b;\
  2436.   register int i,len;\
  2437.   GC_Node;\
  2438. \
  2439.   Check_Type(A,T_Farray);\
  2440.   len = FARRAY(A)->len;\
  2441. \
  2442.   GC_Link(A);\
  2443.   B = farray_make_like(A);\
  2444.   GC_Unlink;\
  2445. \
  2446.   a = FARRAY(A);        /* reassign - A may have moved */\
  2447.   b = FARRAY(B);        \
  2448. \
  2449.   switch(a->type) {\
  2450.   case T_Flonum:\
  2451.     VECLOOP(float,float,inita,initb,statement);\
  2452.     break;\
  2453.   case T_Fixnum:\
  2454.     VECLOOP(int4,int4,inita,initb,statement);\
  2455.     break;\
  2456.   case T_String:\
  2457.     VECLOOP(Zbyte,Zbyte,inita,initb,statement);\
  2458.     break;\
  2459.   default: Panic("v-loopfunc");\
  2460.     break;\
  2461.   }\
  2462.   return B;\
  2463. } /*vloopfunc*/
  2464.  
  2465.  
  2466.  
  2467. /*(DOCENTRY
  2468.   (USAGE "(v-truncate v) => Return integerized vector."
  2469. ))*/
  2470. #define VTRUNCATE    Pvtruncate, "v-truncate", 1,1,EVAL,
  2471.  
  2472. static Object Pvtruncate(A)
  2473.   Object A;
  2474. {
  2475.   Object B;
  2476.   Farray *a,*b;
  2477.   register int i,len;
  2478.   GC_Node;
  2479.  
  2480.   Check_Type(A,T_Farray);
  2481.   len = FARRAY(A)->len;
  2482.  
  2483.   GC_Link(A);
  2484.   B = farray_make(T_Fixnum,len);
  2485.   GC_Unlink;
  2486.   farray_copyshape(A,B);
  2487.  
  2488.   a = FARRAY(A);        /* reassign - A may have moved */
  2489.   b = FARRAY(B);       
  2490.  
  2491.   switch(a->type) {
  2492.   case T_Flonum:
  2493.     VECLOOP(float,int,,, *bp++ = (int)*ap++; );
  2494.     break;
  2495.   case T_Fixnum:
  2496.     VECLOOP(int4,int4,,, *bp++ = (int)*ap++; );
  2497.     break;
  2498.   case T_String:
  2499.     VECLOOP(Zbyte,int,,, *bp++ = (int)*ap++; );
  2500.     break;
  2501.   default: Panic("v-truncate");
  2502.     break;
  2503.   }
  2504.   return B;
  2505. } /*vtruncate*/
  2506.  
  2507.  
  2508.  
  2509. /*(DOCENTRY
  2510.   (USAGE "(v-reverse v) => Return vector in reversed order."
  2511.   "Could be defined as:"
  2512.   "(define (v-reverse a)"
  2513.   "  (let ((size (farray-ref (v-shape a) 0)))"
  2514.   "    (v-[] a (v-- (v-distribute (1- size) size) (v-index size)))))"
  2515. ))*/
  2516. #define VREVERSE    Pvreverse, "v-reverse", 1,1,EVAL,
  2517.  
  2518. LOOPFUNC(Pvreverse,,bp += (len-1);,*bp-- = *ap++;)
  2519.  
  2520.  
  2521. /*(DOCENTRY
  2522.   (USAGE "(v-rotate v a)"
  2523.   "Positive a rotates to the right:"
  2524.   "Result[k+a]:=v[k]."
  2525.   "This could be defined as:"
  2526.   "(define (v-rotate v n)"
  2527.   "  (v-[] v"
  2528.   "    (parlet ((i (v-index (farray-length v))))"
  2529.   "           (v-mod (+ n i) (farray-length v)))))"
  2530. ))*/
  2531. #define VROTATE    Pvrotate, "v-rotate", 2,2,EVAL,
  2532.  
  2533. Object Pvrotate(A,Shift)
  2534.   Object A,Shift;
  2535. {
  2536.   register int4 i,len;
  2537.   int shift;
  2538.   Object B;
  2539.   GC_Node;
  2540.   Error_Tag = "v-rotate";
  2541.  
  2542.   Check_Type(A,T_Farray);
  2543.   len = FARRAY(A)->len;
  2544.   shift = Get_Integer(Shift);
  2545.  
  2546.   GC_Link(A);
  2547.   B = farray_make_like(A);
  2548.   GC_Unlink;
  2549.  
  2550. # define ROTATE(typ_) {\
  2551.   register int4 id = 0 + shift; \
  2552.   register typ_ *ap = (typ_ *)FARRAY(A)->data; \
  2553.   register typ_ *bp = (typ_ *)FARRAY(B)->data; \
  2554.   while (id < 0) id += len;\
  2555.   for( i=0; i < len; i++ ) { \
  2556.     id = id % len; \
  2557.     bp[id] = *ap++; \
  2558.     id++; \
  2559.   }}
  2560.  
  2561.   switch(FARRAY(A)->type) {
  2562.   case T_Fixnum:  ROTATE(int4) break;
  2563.   case T_Flonum:  ROTATE(float) break;
  2564.   case T_String:  ROTATE(char *) break;
  2565.   default: Panic("bad type");
  2566.   } /*switch*/
  2567. # undef ROTATE
  2568.  
  2569.   return B;
  2570. } /*vrotate*/
  2571.  
  2572.  
  2573.  
  2574. /*(DOCENTRY
  2575.   (USAGE "(v-shift v a . fill-value)"
  2576.   "Result[k+a]:=v[k]."
  2577.   "Result has same length as v."
  2578.   "Empty values are filled with fill-value (if provided)"
  2579.   "or with the adjacent array element."
  2580.   "Eventually a boolean fill-value will specify creating a shorter result."
  2581. ))*/
  2582. #define VSHIFT    Pvshift, "v-shift", 0,MANY,VARARGS,
  2583.  
  2584. Object Pvshift(ac,Av)
  2585.   int ac;
  2586.   Object *Av;
  2587. {
  2588.   register int4 i,len,shift;
  2589.   Object B;
  2590.   bool fill;
  2591.   Error_Tag = "v-shift";
  2592.   if ((ac != 2) && (ac != 3)) Primitive_Error("bad # args");
  2593.   Ztrace(("shift nargs=%d\n",ac));
  2594.  
  2595.   Check_Type(Av[0],T_Farray);
  2596.   len = FARRAY(Av[0])->len;
  2597.   shift = Get_Integer(Av[1]);
  2598.   Ztrace(("shift [%d] by %d\n",len,shift));
  2599.  
  2600.   B = farray_make_like(Av[0]);
  2601.  
  2602.   if (ac == 3) {
  2603.     if (TYPE(Av[2]) != FARRAY(Av[0])->type)
  2604.       Primitive_Error("fill type must match array type");
  2605.     fill = TRUE;
  2606.   }
  2607.   else fill = FALSE;
  2608.   Ztrace(("fill = %d\n",fill));
  2609.  
  2610. # define SHIFT(typ_,fill_,value_) {\
  2611.   register typ_ *ap = (typ_ *)FARRAY(Av[0])->data; \
  2612.   register typ_ *bp = (typ_ *)FARRAY(B)->data; \
  2613.   if (shift >= 0) {\
  2614.     register int llen; \
  2615.     if (shift > len) shift = len; \
  2616.     llen = len - shift; \
  2617.     if (fill_)\
  2618.       for( i=0; i < shift; i++ ) bp[i] = value_; \
  2619.     else\
  2620.       for( i=0; i < shift; i++ ) bp[i] = ap[0]; \
  2621.     for( i=0; i < llen; i++ ) { \
  2622.       bp[i+shift] = ap[i]; \
  2623.     }\
  2624.   }\
  2625.   else {\
  2626.     register int llen,len_1; \
  2627.     shift = - shift; \
  2628.     if (shift > len) shift = len; \
  2629.     llen = len - shift; \
  2630.     len_1 = len - 1; \
  2631.     for( i=0; i < llen; i++ )  bp[i] = ap[i+shift]; \
  2632.     if (fill_) \
  2633.       for( i=llen; i < len; i++ ) bp[i] = value_; \
  2634.     else\
  2635.       for( i=llen; i < len; i++ ) bp[i] = ap[len_1]; \
  2636.   }}
  2637.  
  2638.  
  2639.   switch(FARRAY(Av[0])->type) {
  2640.   case T_Fixnum: {
  2641.     int value;
  2642.     if (fill) value = Get_Integer(Av[2]);
  2643.     SHIFT(int4,fill,value)
  2644.     break;
  2645.   }
  2646.   case T_Flonum: {
  2647.     float value;
  2648.     if (fill) {
  2649.       Check_Type(Av[2],T_Flonum); value = FLONUM(Av[2])->val;
  2650.     }
  2651.     SHIFT(float,fill,value)
  2652.     break;
  2653.   }
  2654.   case T_String: {
  2655.     int value;
  2656.     if (fill) {
  2657.       Check_Type(Av[2],T_Character);      value = CHAR(Av[2]);
  2658.     }
  2659.     SHIFT(char,fill,value) 
  2660.     break;
  2661.   }
  2662.   default: Panic("bad type");
  2663.   } /*switch*/
  2664. # undef SHIFT
  2665.  
  2666.   return B;
  2667. } /*vshift*/
  2668.  
  2669.  
  2670.  
  2671. /*(DOCENTRY
  2672.   (USAGE "(v-head v d)"
  2673.   "Return the vector resulting from discarding the d lowest elements"
  2674.   "of v."
  2675.   "(v-head (% 1 2 3) 1) => (% 2 3)"
  2676. ))*/
  2677. #define VHEAD    Pvhead, "v-head", 2,2,EVAL,
  2678.  
  2679. Object Pvhead(A,Discard)
  2680.   Object A,Discard;
  2681. {
  2682.   int discard,len,hlen;
  2683.   Object B;
  2684.   GC_Node;
  2685.   Error_Tag = "v-head";
  2686.  
  2687.   Check_Type(A,T_Farray);
  2688.   discard = Get_Integer(Discard);
  2689.   len = FARRAY(A)->len;
  2690.   if (discard >= len)   Primitive_Error("array is not that long");
  2691.   if (discard < 0) Primitive_Error("bad arg (negative)");
  2692.  
  2693.   hlen = len - discard;
  2694.  
  2695.   GC_Link(A);
  2696.   B = farray_make(FARRAY(A)->type,hlen);
  2697.   GC_Unlink;
  2698.  
  2699.   switch(FARRAY(A)->type) {
  2700.   case T_Flonum:
  2701.     Zbcopy( (char *) (((float *)FARRAY(A)->data)+discard),
  2702.            (char *)FARRAY(B)->data, hlen * sizeof(float));
  2703.     break;
  2704.   case T_Fixnum:
  2705.     Zbcopy( (char *) (((int *)FARRAY(A)->data)+discard),
  2706.            FARRAY(B)->data, hlen * sizeof(int4));
  2707.     break;
  2708.   case T_String:
  2709.     Zbcopy( ((char *)FARRAY(A)->data)+discard,
  2710.            FARRAY(B)->data, hlen * sizeof(char));
  2711.     break;
  2712.   default: Panic("v-head");
  2713.   }
  2714.  
  2715.   return B;
  2716. } /*head*/
  2717.  
  2718.  
  2719.  
  2720. /*(DOCENTRY
  2721.   (USAGE "(v-tail v d)"
  2722.   "Return the vector resulting from discarding the d highest elements"
  2723.   "of v."
  2724.   "(v-tail (% 1 2 3) 1) => (% 1 2)"
  2725. ))*/
  2726. #define VTAIL    Pvtail, "v-tail", 2,2,EVAL,
  2727.  
  2728. Object Pvtail(A,Discard)
  2729.   Object A,Discard;
  2730. {
  2731.   int discard,len,hlen;
  2732.   Object B;
  2733.   GC_Node;
  2734.   Error_Tag = "v-tail";
  2735.  
  2736.   Check_Type(A,T_Farray);
  2737.   discard = Get_Integer(Discard);
  2738.   len = FARRAY(A)->len;
  2739.   if (discard >= len) Primitive_Error("array is not that long");
  2740.   if (discard < 0) Primitive_Error("bad arg (negative)");
  2741.  
  2742.   GC_Link(A);
  2743.   B = farray_make(FARRAY(A)->type,len-discard);
  2744.   GC_Unlink;
  2745.  
  2746.   hlen = len - discard;
  2747.  
  2748.   switch(FARRAY(A)->type) {
  2749.   case T_Flonum:
  2750.     Zbcopy( (char *)FARRAY(A)->data, FARRAY(B)->data, hlen * sizeof(float));
  2751.     break;
  2752.   case T_Fixnum:
  2753.     Zbcopy( (char *)FARRAY(A)->data, FARRAY(B)->data, hlen * sizeof(int4));
  2754.     break;
  2755.   case T_String:
  2756.     Zbcopy( (char *)FARRAY(A)->data, FARRAY(B)->data, hlen * sizeof(char));
  2757.     break;
  2758.   default: Panic("v-tail");
  2759.   }
  2760.  
  2761.   return B;
  2762. } /*tail*/
  2763.  
  2764.  
  2765.  
  2766. /*(DOCFINIT)
  2767. */
  2768. /*%%%%%%%%%%%%%%%% link %%%%%%%%%%%%%%%%*/
  2769.  
  2770. static struct primdef Prims[] = {
  2771. /* multi-dimensional arrays */
  2772.   VARRAY
  2773.   VARRAYSET
  2774.   VARRAYREF
  2775.  
  2776. /* shape/reference functions */
  2777.   VSHAPE
  2778.   VRAVEL
  2779.   VRESHAPE
  2780.   VTRANSPOSE
  2781.   VCOMPRESS
  2782.   VAPPEND
  2783.   VSELECT
  2784.  
  2785.   VMAPCOUNT
  2786.   VREFERENCE    /* cf gather/scatter */
  2787.   VREFERENCEb
  2788.  
  2789.   VGATHER
  2790.   VSCATTER
  2791.  
  2792. /* monadic functions */
  2793.  
  2794.   VSIN
  2795.   VCOS
  2796.   VSQRT
  2797.   VEXP
  2798.  
  2799.   VABS
  2800.   VNOT
  2801.  
  2802.   VRANDOM
  2803.  
  2804. /* return highest type */
  2805.   VMUL
  2806.   VADD
  2807.   VSUB
  2808.   VDIV
  2809.   VMIN
  2810.   VMAX
  2811.   VMOD
  2812.   VIMOD /*mod as %, returns integer */
  2813.  
  2814. /* return bool vector */
  2815.   VEQ
  2816.   VNE
  2817.   VLE
  2818.   VLT
  2819.   VGT
  2820.   VGE
  2821.   VAND
  2822.   VOR
  2823.   VXOR
  2824.  
  2825. /* return float regardless */
  2826.   VPOW
  2827.   VATAN2
  2828.  
  2829. /* return scalar same type as a[0] */
  2830.   VPLUSREDUCE
  2831.   VPLUSREDUCEb
  2832.   VPLUSREDUCEc
  2833.   VMINUSREDUCE
  2834.   VMINUSREDUCEb
  2835.   VMINUSREDUCEc
  2836.   VMULREDUCE
  2837.   VMULREDUCEb
  2838.   VMULREDUCEc
  2839.   VMINREDUCE
  2840.   VMINREDUCEb
  2841.   VMINREDUCEc
  2842.   VMAXREDUCE
  2843.   VMAXREDUCEb
  2844.   VMAXREDUCEc
  2845.  
  2846. /* return same type of vector */
  2847.   VPLUSSCAN
  2848.   VPLUSSCANc
  2849.   VMINUSSCAN
  2850.   VMINUSSCANc
  2851.   VMULSCAN
  2852.   VMULSCANc
  2853.   VMINSCAN
  2854.   VMINSCANc
  2855.   VMAXSCAN
  2856.   VMAXSCANc
  2857.  
  2858. /* misc */
  2859.   VDISTRIBUTE
  2860.   VINDEX
  2861.   VTRUNCATE
  2862.   VREVERSE
  2863.   VROTATE
  2864.   VSHIFT
  2865.   VHEAD
  2866.   VTAIL
  2867.  
  2868.   VTEST
  2869.  
  2870.   (Object (*)())0, (char *)0, 0,0,EVAL
  2871. };
  2872.  
  2873.  
  2874. void Init_vector()
  2875. {
  2876.   ZLprimdeftab(Prims);
  2877. }
  2878.  
  2879. #endif /*ELKVECTOR*/
  2880.